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CHAPTER  I 


INTRODUCTION 


Many  important  problems  of  physical  geodesy  are  being  solved  by 
integrals  extended  over  the  entire  Earth.  Some  examples  are  Stokes’ 
and  Vening-Meinesz*  formulae  [Heiskanen  and  Moritz,  1967]  and 
Molodenskii’s  solution  of  the  boundary  value  problem  of  physical 
geodesy  [Moritz,  1980].  The  solution  of  these  problems  is  formulated 
assuming  that  gravity  is  known  everywhere  on  the  Earth’s  surface. 
However,  this  is  hardly  ever  the  case.  Therefore,  there  is  a  well 
justified  need  for  gravity  interpolation  and  extrapolation. 

The  scope  of  this  study  is  to  investigate  two  deterministic  methods 
for  gravity  field  approximation.  The  predictors  are  deterministic  in  the 
sense  that  no  stochastic  properties  of  the  gravity  field  are  involved. 
On  the  other  hand,  one  feature  of  both  methods  is  that  any  linear 
functional  of  the  disturbing  potential  can  be  used  as  observable  and/or 
quantity  to  be  predicted. 

The  first  method  was  proposed  by  Bjerhammar  [1964].  He  assumed 
that  observations  are  given  at  a  finite  number  of  stations.  The 
disturbing  potential  is  assumed  harmonic  (i.e.,  it  satisfies  Laplace’s 
equation)  down  to  a  sphere  fully  internal  to  the  Earth.  The  gravity 
anomalies  at  the  surface  of  the  internal  sphere  are  solved  for  by  a 
downward  continuation  process  and  then  they  are  used  to  perform 
predictions  by  an  upward  continuation  integral.  This  formulation  and 
solution  of  the  discrete  geodetic  boundary  value  problem  is  treated  in 
Chapter  2. 

Secondly,  Hardy’s  biharmonic  potential  technique  is  considered. 
According  to  this  approach  the  disturbing  potential  can  be  shown  to 
satisfy  the  biharmonic  equation.  The  biharmonic  sources  that  generate 
the  disturbing  potential  can  be  estimated  based  on  observations  and 
then  they  can  be  used  to  predict  gravity  related  quantities.  The 
derivation  of  the  alternate  integral  for  the  disturbing  potential,  the 
biharmonic  equation  and  its  solutions  in  terms  of  spherical  harmonics, 
as  well  as  expressions  of  the  gravity  anomalies  and  the  deflections  of 
the  vertical  in  terms  of  the  biharmonic  sources  are  given  in  Chapter  3. 

Chapter  4  includes  a  detailed  description  of  the  terrain  and  the 
data  coverage  as  well  as  data  reduction  computations  for  the  White 
Sands  Test  Area  in  New  Mexico.  The  White  Sands  test  data  were  used 
to  test  both  predictors. 
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Thia  work  is  continued  with  a  detailed  account  of  the  teats 
performed  with  both  methods.  Variations  of  the  methods  are  tested 
and  the  results  are  discussed  and  summarized.  A  comparison  is 
attempted  with  the  results  of  the  four  methods  tested  at  the  same  area 
[Kearsley  et  al.,  1985].  The  aforementioned  testa  are  given  in  Chapter 
5. 

Finally ,  a  summary  of  this  investigation  together  with  the 
conclusions  drawn  and  recommendations  for  future  work  is  given  in 
Chapter  6. 


CHAPTER  II 


THE  BJERHAMMAR  PROBLEM 


2.1  Introduction 


An  explicit  solution  to  the  geodetic  boundary  value  problem  was 
published  by  George  Gabriel  Stokes  in  1849  [Heiskanen  and  Moritz, 
1967,  p.  94].  The  underlying  assumptions  of  his  solution  are  that  the 
mathematical  figure  of  the  Earth  (the  geoid)  is  approximated  by  a 
sphere;  there  are  no  masses  outside  the  geoid;  and  that  gravity  is 
known  everywhere  on  the  geoid.  However,  gravity  is  measured  on  the 
surface  of  the  Earth,  and  the  application  of  Stokes'  formula  requires 
on  one  hand  the  reduction  of  measured  gravity  down  to  the  geoid  and 
on  the  other  hand  the  absense  of  masses  above  the  geoid  [ibid,  p. 
126].  In  order  to  solve  these  two  problems  gravity  reductions  must  be 
used.  These  gravity  reductions  not  only  require  an  assumption  for  the 
density  of  the  external  masses  [Bjerhammar,  1964,  p.  9],  they  also 
introduce  a  change  of  the  geoid,  the  indirect  effect  [Heiskanen  and 
Moritz,  1967,  p.  141].  Molodenskii  in  1945  stated  that  the  geoid  cannot 
be  determined  without  knowledge  of  the  mass  distribution  outside  it 
[Bjerhammar,  1963,  p.  3]. 


Malkin  in  1939  redefined  the  principal  problem  of  physical  geodesy 
so  that  the  physical  surface  of  the  Earth  became  the  unknown  [ibid,  p. 
9]. 


In  1948,  Molodenskii  presented  the  solution  to  the  problem  of 
determining  the  physical  surface  of  the  Earth  from  gravity 
measurements  [Moritz,  1980,  p.  330],  [Bjerhammar,  1964,  p.  3].  The 
problem  is  non-linear  and  the  notion  of  the  telluroid  is  introduced  for 
its  Taylor  linearisation  [Moritz,  1980,  p.  337].  Molodenskii’s  solution, 
Brovar’s  solution  and  a  solution  by  analytical  continuation  are  given  by 
Moritz,  [1980,  pp.  354-388].  Bjerhammar  [1963,  p.  7]  gave  credit  to 
Molodenskii  for  his  elegant  solution  for  the  distrubing  potential  but  he 
also  realized  that  the  problem  of  Molodenskii  is  continuous 
[Bjerhammar,  1975,  p.  185],  [Bjerhammar,  1964,  p.  3]  and  the 
observations  are  only  made  at  discrete  points. 

As  a  result,  Bjerhammar  defined  the  discrete  geodetic  boundary 
value  problem  as  follows  [Bjerhammar,  1975,  p.  185],  [Bjerhammar,  1964, 
p.  14],  [Bjerhammar,  1963,  p.  17],  [Moritz,  1980,  p.  95],  [Bjerhammar, 
1986,  p.  1]:  A  finite  number  of  gravity  data  is  given  for  a 
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non-spherical  surface  and  it  is  required  to  find  such  a  solution 
(gravity  field)  that  the  boundary  values  for  the  gravity  data  are 
satisfied  at  all  given  points.  The  discrete  boundary  value  problem 
avoids  the  singularities  of  downward  continuation  [Moritz,  1980,  p.  95] 
and  has  no  need  for  uniform  approximation  theorems  (such  as  Runge’s 
or  Keldych-Lavrentieff’s)  [Bjerhammar,  1973,  p.  480].  Stokes  solved  the 
geodetic  boundary  value  problem  assuming  continuous  data  coverage  on 
the  boundary.  An  error  is  introduced  in  the  actual  implementation  of 
his  solution  because  an  integral  is  being  approximated  by  a  finite  sum. 
On  the  other  hand,  Bjerhammar  realized  the  fact  that  there  is  only  a 
finite  amount  of  data  and  consequently  he  formulated  the  discrete 
geodetic  boundary  value  problem. 


2.2  Formulation  of  the  Problem 
2.2.1  Definitions 

This  study  will  follow  well  established  geodetic  practices  in  which 
the  geodetic  boundary  value  problem  is  independent  of  time  and  the 
space  outside  the  boundary  is  empty  [Moritz,  1980,  p.  330].  This  is  to 
say  that  the  atmospheric  and  the  tidal  effects  have  been  taken  into 
account  by  corrections  to  the  observed  gravity  [ibid,  pp.  425,  330], 
[Sanao’,  1981,  p.  13].  The  Earth  is  assumed  to  be  a  rigid  body  which 
rotates  with  constant  angular  velocity  around  a  fixed  axis  which  passes 
through  its  center  of  mass  [Moritz,  1980,  pp.  330,  447],  [Sanso’,  1981,  p. 
13].  An  Earth-fixed  rectangular  Cartesian  coordinate  system  is  defined 
such  that  the  origin  is  at  the  center  of  mass  of  the  Earth.  Its  z  axis 
coincides  with  the  Earth’s  mean  axis  of  rotation,  the  xz  plane  coincides 
with  the  mean  Greenwich  meridian  plane  and  the  y  axis  is 
perpendicular  to  the  xz  plane  such  that  the  system  is  right-handed 
[Moritz,  1980,  p.  2J.  Therefore,  the  figure  and  gravity  field  of  the 
Earth  as  well  as  the  coordinate  Bystem  are  assumed  to  be  constant  in 
time  [ibid,  p.  477]. 

Let  W  be  the  actual  gravity  potential  of  the  Earth  and  let  U  be  a 
norsial  gravity  potential  which  is  an  analytic  approximation  to  W;  U  is 
usually  taken  as  the  potential  of  an  equipotential  ellipsoid  [Moritz, 
1980,  p.  337].  Let  us  also  denote  by  2  the  actual  gravity  vector  and 
by  ?  the  normal  gravity  vector.  They  are  defined  as  follows  [ibid,  p. 
337]: 


2  -  gradW, 

(2-1) 

f  -  gradU. 

(2-2) 

A  point  P  of  the  geoid  is  projected  onto  the  point  Q  of  the 
ellipsoid  by  means  of  the  ellipsoidal  normal  (Figure  1).  The  gravity 
anomaly  vector  &2  is  defined  as  the  difference  [Heiskanen  and  Moritz, 
1967,  p.  83], 
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=  ip  "  tn 

where  $p  is  the  actual  gravity  vector  at  P  and  the  normal  gravity 
vector  at  Q.  The  difference  in  magnitude 


Ag  =  gp  -  ya 

is  the  gravity  anomaly  [ibid,  p.  83). 


(2-3) 


«**U  (WaW.) 


rtftrtacc  ,„_w, 


Figure  1.  Definition  of  the  gravity  anomaly. 


The  disturbing  or  anomalous  potential  T  is  defined 
12] 

T  =  W  -  U. 


[Moritz,  1980,  p. 


(2-4) 


The  fundaaental  equation  of  physical  geodesy,  neglecting 
i2  _  ifi  +  0(e4) 

which  is  smaUer  than  0.5  mgal  [Cruz,  1985],  is  [Heiskanen  and  Moritz, 
1967,  pp.  86,  91] 


7S  -  ^  -«-«*• 


(2-5) 


Introducing  the  usual  spherical  approximation,  i.e.,  the  tolerance  of 
a  relative  error  in  the  order  of  3x10" 3  in  equations  relating  quantities 
of  the  anomalous  field,  one  can  write  [Heiskanen  and  Morits,  1967,  pp. 
87-88],  [Bjerhammar,  1986,  p.  5]: 


.  <M  .  » 

7  ra  1  Jh 

and  therefore 


_  Lii  =  _  A.ix  =  _ 

7  Jh  y  Jr  y  JrlraJ 
Substituting  in  (2-5)  one  gets 


JT  2T 
Jr  r 


+ 


Ag 


0. 


(2-6) 


2  .  3T  JT 

7  and  ah  =  a7 


(2-7) 


2.2.2  Pizzetti’a  Formula 

Let  us  now  introduce  a  sphere  <r  fully  internal  to  the  Barth,  of 
radius  r0  and  with  its  center  at  the  center  of  mass  of  the  Earth.  This 
sphere  is  hereafter  called  the  Bjerhammar  sphere,  the  internal  sphere, 
or  the  geosphere. 

The  aim  is  to  solve  the  discrete  boundary  value  problem,  i.e.,  given 
a  finite  number  of  observed  gravity  data,  to  find  a  solution  (disturbing 
potential)  with  the  following  properties  [Bjerhammar,  1986,  p.  8]: 

(a)  T  is  harmonic  outside  <r, 

(b)  T  is  regular  at  infinity,  i.e.,  ^b(Tt)  =  constant,  and 

(c)  All  observations  are  satisfied. 

If  one  aaaumes  no  masses  outside  the  sphere  a  and  denote  the 
gravity  anomalies  on  <r  by  Ag*  then  the  disturbing  potential  T  is  given 
by  Pissetti's  formula  [Heiskanen  and  Moritz,  1967,  p.  93],  [Bjerhammar, 
1986,  p.  7]: 

T(r,  ♦,  X)  =  JJs(r,  «)Ag*de  (2-8) 

where  [Heiskanen  and  Moritz,  1967,  p.  93] 


S(r,  *)  •  ^  -  &  cm.  (e  ♦  3<n 

with  cornu  =  sin+sintf  +  cos4cos*,cos(X-X( ) 
and  i*  =  ra  +  r0*  -  2r0rcos«; 


r-r„coa«,+<) 

2r  J 


(2-9) 

(2-10) 

(2-11) 
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also,  ♦,  X  are  the  geodetic  latitude  and  longitude  respectively  and  r  is 
the  geocentric  distance  of  the  point  at  which  T  is  computed  and  ♦, ,  X, 
are  the  geodetic  latitude  and  longitude  respectively  of  the  infinitesimal 
surface  element  d<r. 

Let  £  =  [♦  A]r.  Let  us  also  define  the  evaluation  [Moritz,  1980,  pp. 
37,  200]  or  Dirac  "delta”  functional  6  [Bjerhammar,  1986,  p.  14]  by  the 
following  relation  (*  is  the  geosphere,  R  is  the  set  of  real  numbers  and 
f  is  some  function  on  <r): 

f  *  *  -*♦  H:  ^  JJ  f (£)«(£  -  x,)de  =  f(x,)  (2-12) 

The  Ag*  in  equation  (2-8)  is  rewritten  as  [Bjerhammar,  1986,  p.  14]: 


Ag*(x)  =  I  Ag^x,  )5(x-x,)  (2-13) 

i=» 

The  Ag*(£| )  values  are  a  set  of  fictitious  anomalies  that  generate  the 
observations  at  the  given  points  [Sjoberg,  1978,  p.2],  [Katsambalos, 
1981,  p.58].  The  basic  postulate  of  the  Dirac  Impulse  method  is  that 
Ag*(x)  =  0  everywhere  on  the  geosphere  with  the  exception  of  n  points 
associated  to  the  given  observations  [Bjerhammar,  1986,  p.14]. 

Equation  (2-13)  is  the  mathematical  formulation  of  the  basic  postulate. 

Substituting  (2-13)  in  (2-8)  one  obtains 


T(r,  *)  =  JJs(r,  x,  x„)  Ag*(x„)de„  = 

a 

=  JJ  S(r,  x,  xN)  •^Ag*(x1)A(xs-x1)dffM  = 
=  rot|tAg*(x,)-  ^  JJs(r,  x,  xM)«(xH-x, )de„ 

If 


which  by  (2-12)  yields  (the  subscript  M  in  the  above  derivation  was 
merely  introduced  to  denote  the  moving  point  M  on  the  sphere): 


T(r,4,X)  =  r0 


X»)-S(r,4, A, ♦,,>,) 


Introducing  t  and  d  by 


(2-14) 


(2-15) 


8 


d  =  {1  +  ta  -  2tcos«}*  (2-16) 

one  gets  from  (2-11) 

t  =  rd  (2-17) 

Letting 

u  =  |  (1  -  tcosu  +  d)  (2-18) 

the  generalized  Stokes'  function  S(r,  a)  in  (2-9)  becomes 

S(r,  »)  =  t(l  -  3d  +  j  -  t cos« (5  +  31n  u))  (2-19) 

Therefore,  (2-14)  becoaes 

T(r,t,X)  =  r0t  2  (1  -  3d  +  \  -  tcosw(5  +  3inu))Agf  (2-20) 

i=i  a 


2.2.3  Gravity  anomalies  and  vertical  deflections  in  terms  of  Dirac 
anomalies 

All  the  quantities  related  to  the  disturbing  potential,  i.e.,  all  its 
linear  functionals  auch  as  gravity  anomalies  or  vertical  deflections  can 
be  evaluated  by  applying  the  pertinent  linear  operators  to  T  in  (2-20). 

(a)  Gravity  Anomalies: 

From  (2-7)  one  sees  that  Ag  =  LT,  where 


AT 

In  order  to  compute  yp  ,  the  following  derivatives  are  needed 

Ad  _  t-coso 
At  "  d 


Au  If t-coso 
At  "  2l  d 


CO»(J 


(2-21) 

(2-22) 


as  T(r,  ♦,  X)  =  r0  |  b,Agf,  with 

f =i 


Rewriting  (2-20) 
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b«  =  t  -  3dt  +  j - 5tacos6>  -  3t JcoB<ainu 


(2-23) 


one  has 


A«f. 


(2-24) 


Now, 


3b,  9t 
at  ar  • 


and 


at 

ar 


Thus, 

—  .UP  g  **-?,!.  ».*■ 
-  '  r«  ,1.  IP  (-  ?M  -  2t  ,1.  «>•«* 


(2-25) 


or, 


Ag  =  ~ta  2  M,Agf  (2-27) 

i  =  » 

2b,  ab, 

where  M,  =  -^  - 

The  M(  coefficients  can  be  further  reduced  (see  Appendix  A.1)  to  yield 
Mt  =  1  +  “jpP'  +  3tcos« 
so  that  (2-27)  becoaes 

kg  -  -  t*  |Jl  +  +  3tcos«) Agf 

or,  rearranging  terns 

Ag  -,1,^  ~  3t*co^  ~  ta)A*f  (2-28) 

which  is  of  course  the  same  equation  Bjerhammar  came  up  with  by 
discretizing  Poisson’s  equation  (Bjerhammar,  1986,  p.8], (Bjerhammar, 
1964,  p.20]. 

(b)  Deflections  of  the  vertical: 

The  vertical  deflections  north  {{)  and  east  (ij)  are  given  by 
Heiakanen  and  Moritz  [1967,  p.235]: 
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*  _  1  .  .  _  1  - 
(  =  -  -  V  ~  -  -  «x. 

where  [ibid,  p.233] 

x  -  1  «T  .  .  _  1  £T 
*♦  '  r  if  *  6x  “  rcos*  ax 

or,  upon  substitution 


,  i_  rr 

*  '  yr  if 


,  _  _  _ _  3T 

l1*  ~  7rcost  i\ 

On  the  other  hand 


i2  i"  ,  3T  3T  3u 
if  =  iu  if  and  ax  =  iu  ax 


(2-29) 


and  [ibid,  p.234] 


iu 


=  -cosd  and  ^  =  -cos ♦sin* 

therefore 

ft]  ,  fcose]  1_  IT 
lifi  vsinoti  yr  iu 


(2-30) 


Mow  -  r0t  S  A|Xgf,  with  At  =  [l  -  3d  +  ^  -  5tcos«-3tcos«Anu] 
The  A(  derivatives  can  be  evaluated  (see  Appendix  A. 2)  to  yield 
A,  =  tain*  (8  -  +3Jnu] 


thus 

£  ■  r.fj.js  -  \s  -  *  3<m.).inaigf 

Substitution  of  (2-31)  in  (2-30)  yields 

{*}-  - 1*  -  ^  (-M 

Finally,  using  Heiskanen  and  Noritz  [1967,  p.113] 


(2-31) 


(2-32) 
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(sinwcoset  =  cos+sin+i  -  s in+cost < cos ( X ^ -X ) ) 

sinwsinci  =  cos*,sin(A,-A)  i  x  ' 

one  gets 

m  ■  r,U8  -*-3-w-+  *h  e^£ff“*,~(x,-wM 

(2-34) 

Equations  (2-28)  and  (2-34)  can  be  used  both  to  compute  Agf  from 
observations  and  to  predict  Ag,  (,  ij. 


2.2.4  A  truncated  sum 

Given  a  continuous  function  f(x)  and  a  constant  x0,  one  can  form  a 
sequence  fn  =  f(nx0),  n  «  N.  This  operation  is  called  sampling  and  Xo 
is  called  the  sampling  interval.  Conversely,  given  a  sequence  fn,  one 
can  construct  a  continuous  function  f(x),  by  an  operation  called 
pulsing  and  defined  as 


f(x)  =  I  f„0(x  -  nx0)  =  I  f(nx#)«(x  -  nx0) 

n=0  n=o 

where  6  is  the  Dirac  delta  functional.  The  aforementioned  operations 
are  well  established  in  the  analysis  of  Linear  Systems  and  can  be 
found  in  many  electrical  or  mechanical  engineering  texts  such  as  Brown 
[1961,  p.  176],  Aseltine  [1958,  p.  247]  and  Tretter  [1976,  p.  85]. 

Bquation  (2-13)  is  the  two-dimensional  analog  of  the  above  equation 
defining  pulsing.  The  only  difference  is  that  pulsing  as  defined  above 
is  the  superposition  of  an  infinite  number  of  impulses  whereas  (2-13) 
represents  a  finite  sum. 

Bjerhammar  realised  that  in  applications  one  has  only  a  finite 
number  of  data,  and  therefore  one  can  only  solve  for  a  finite  number 
of  impluses.  Consequently  he  truncated  the  sum  in  (2-13)  at  some 
integer  number  n  equal  to  the  number  of  observations. 

However,  only  an  infinite  number  of  Ag*  values  is  able  to 
recontruct  the  original  signal  (recall  that  in  order  to  recover  T  one 
needs  an  infinite  number  of  harmonics).  Therefore,  the  gravity  related 
quantities  predicted  from  n  Ag*  values  will  naturally  not  include  the 
contribution  of  the  truncated  terms  beyond  n.  This  is  to  say  that 
equations  such  as  (2-20),  (2-28)  and  (2-34)  can  only  be  considered  as 
approximations. 
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2.3  Symmetric  Kernel  Approach 

In  the  original  approach  the  Dirac  Impulses  were  located  on  the 
surface  of  the  sphere  a .  In  the  symmetric  kernel  approach  the  Dirac 
Impulses  are  located  at  positions  (Figure  2)  with  geocentric  radii 
inversely  proportional  to  the  geocentric  radii  of  the  observations,  i.e. 
the  Dirac  Impulse  corresponding  to  the  observation  at  r<  is  located  at 
r0«,  where  [Bjerhammar,  1986,  p.  48] 

(2-35) 


Figure  2.  Symmetric  Kernel  Approach. 


Bvery  observation  will  be  associated  with  a  Dirac  Impulse  at  r0) . 
Considering  n  observations  one  will  have  to  consider  n  Dirac  Impulses 
located,  in  principal,  at  n  different  geospheres.  The  disturbing 
potential  associated  with  the  i-th  Impulse  will  be  given  by  (see 
equation  (2-14)) 


T,(r,  ♦,  X)  =  r0,S(r,  ♦,  X,  X,)Agf  (2-36) 

and  ATj  =  0  at  r  >  r0t.  The  total  potential  will  be 
T(r,  ♦,  X)  -  ^ro)S(r,  ♦,  X,  4j,  X , )  A  gf 


(2-37) 
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by  the  superposition  principle  of  the  potential  [Heiskanen  and  Moritz, 
1967,  p.2) 


and  AT  =  0  at  r  >  max(r01,  _ _  r0„) 


min(rt,  ....  rn) 

Introducing  t  by  t  =  and  d  by  (2-16)  one  has 

t  =  lax  =  Eal 
r  rr  f 

and 


(2-38) 


S(r,  «)  =  t(l  -  3d  +  5  -  tcosu(5  +  3inu)) 

a 


(2-39) 


Therefore  (2-37)  becones 

T(r,4,X)  =  £  r0tt(l  -  3d  +  \  -  tcosu(5  +  3*nu))Agf  (2-40) 

i  =  i  a 

Similarly  as  in  the  original  approach  one  derives 

Ag  =  l  _  3t3cos«  -  t*  j Agf  (2-41) 

{ ?  ]  a  1 13(8._  2  _  3(±ni  +3inu]{cosJsin^:si"tcos*,cos(x,_x)l4gf 

1  ij  >  y  ,  =  i  1  d3  2ud  Hcos4,sin(A,-A)  J 

(2-42) 


2.4  The  Linear  System 

If  one  assumes  n  observations  comprising  the  vector  I  then  the 
linear  system  is 

i  =  GAg*  (2-43) 

In  general  the  vector  of  the  Dirac  Impulses  Ag*  will  be  of 
dimension  m.  The  elements  of  G  can  be  taken  from  (2-28)  or  (2-34)  for 
or  (,  rj  observations  in  the  original  approach  and  from  (2-41),  (2-42) 
for  the  same  type  of  observations  in  the  symmetric  kernel  approach. 
Note  that  G  is  of  full  row  rank  as  long  as  there  are  no  two 
observations  l<,  I  j  of  the  same  type  (e.g.,  Ag  or  ()  with  r(  =  rj. 

There  are  three  distinct  possibilities: 
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(i)  Exact  solution  (m  =  n): 

The  Dirac  Impulses  and  their  covariance  matrix  respectively  are 
given  by 

Ag*  =  <S~li  (2-44) 


Ug*  =  <r‘Ij<rT  (2-45) 

(ii)  Overdetermined  case  (m  <  n): 

In  this  case  the  traditional  least-squares  method 
applies  and  the  solution  is 

Ag*  =  (GTEi1G)~16TIj*  4 

Ua*  =  WU'G)-' 

The  a  posteriori  variance  of  unit  weight  a%  is  given  by 

H  ,  ... 

n  -  » 

where  <r£  is  the  a  priori  variance  of  unit  weight. 

It  should  be  kept  in  mind  that  the  property  of  the  reproducibility 
of  the  observations  is  lost  in  this  case. 

1  Underdetermined  case  (m  >  n): 

In  this  case  equation  (2-43)  represents  a  system  of  equations  with 
an  infinite  number  of  solutions  since  G  has  full  row  rank  [Dermanis, 
1985,  p.  215].  From  this  infinite  number  of  solutions,  one  particular 
solution  may  be  chosen  by  the  minimum-norm  criterion  Ag*T-Ag*  =  min. 
This  solution  is  unique  in  the  sense  that  the  norm  of  any  other 
solution  is  at  least  as  large  as  the  norm  of  Ag*  or  larger.  In  order  to 
find  the  solution  that  minimizes  the  norm  of  Ag*  and  satisfies  (2-43) 
one  needs  to  find  Ag*  that  satisfies 

♦  =  Ag*TAg*  +  2kT(GAg*  -  4)  =  ain  (2-49) 

The  ainiaua  of  4  with  respect  to  Ag*  is  attained  at 

Ujg  =  0  <=>  2Ag*T  +  2kTG  =  0  <=>  Ag*  =  GTk  (2-50) 

On  the  other  hand 


[Uotila,  1985) 

(2-46) 

(2-47) 

(2-48) 


0  =  GAg*  -  J  =  GGTk  -  4  <=>  k  =  (GGT)-‘4,  thus 


15 


Ag*  =  GT(GGT)-‘A  (2-51) 

Obviously  G  is  (nxm)  with  rank  n  (full  row  rank).  Now  GGT  is 
(nxn)  with  the  same  rank  as  G,  thus  GGT  is  of  full  rank  and  therefore 
invertible.  The  covariance  matrix  of  Ag*  can  be  computed  from  the  law 
of  propagation  of  covariances  (Uotila,  1985,  p.  4]: 

£a9*  =  G^GG*)"*  Zl  (GT(GGT)_1  )T  =  GT(GGT)-*  it  (GGT)-1G,  thus 
lAg*  =  GT(GGT)_1  Zj  (GGT)-‘G  (2-52) 


It  is  trivial  to  notice  that  the  solution  (2-51)  is  the  pseudo  inverse 
solution.  For  a  full  row  rank  matrix  G  one  has  (Uotila,  1982b] 

G+  =  Gt(QGt)-1  '  (2-53) 

thus  (2-51)  and  (2-52)  can  be  written  as 

Ag*  =  G+J  (2-54) 

Ia9*  =  G4  Zl  (G+)t  (2-55) 

The  uniqueness  of  the  pseudo  inverse  quarantees  that  Ag*  in  (2-54)  is 
unique. 


2.5  Propagation  of  Data  Noise  into  the  Predictions 

Let  us  assume  that  we  will  perform  predictions  at  q  stations.  The 
vector  f  of  predictions  will  be  at  dimension  q  and  it  will  be  given  by 

f  =  RAg*  (2-56) 

where  the  elements  of  the  (qxm)  matrix  R  will  be  given  by  (2-28)  for 
Ag  predictions  and  by  (2-34)  for  ((,*))  predictions.  The  covariance 
matrix  E;  of  the  prediction  vector  f  will  be  given  by  [Uotila,  1985,  p.4] 

Ef  =  R‘Slgr*T  (2-57) 

where  E&g*  is  the  covariance  matrix  of  Ag*  as  computed  from  (2-45)  for 
the  exact  case,  (2-47)  for  the  overdetermined  case  and  (2-55)  for  the 
underdetermined  case. 

2.6  On  the  Location  of  the  Dirac  Impulses 

From  the  theoretical  standpoint  there  is  no  reason  to  prefer  any 
location  over  any  other  for  the  location  of  the  Dirac  Impulses.  From 
the  numerical  standpoint  however,  one  should  prefer  the  location  that 
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yields  optimal  stability  in  the  sense  of  maximal  diagonal  dominance  of 
the  matrix  to  be  inverted  as  done  in  Bjerhammar  [1986,  p.  30].  The 
matrix  G  in  (2-44)  will  be  selected  to  investigate  the  location  issue 
since  this  issue  is  only  numerically  and  not  theoretically  relevant.  The 
resulting  location  is  going  to  be  used  regardless  of  the  observation 
type  and  of  the  solution  type  (i.e.»  exact  solution,  adjustment  solution 
or  minimum  norm  solution).  The  development  will  follow  Bjerhammar's 
ideas  [ibid,  p.  30]. 

Let  the  Dirac  Impulse  Ag*  corresponding  to  some  observation  Ag  at 
r,  4,  X  be  located  at  r0,  40>  X0. 


The  corresponding  element  of  G  in  the  main  diagonal  will  be  (see 
eq.  (2-28)) 


g  =  -  3t3coso  -  ta  (2-58) 
with  d  =  {1  +  t2  -  2tcos (2-59) 
and  cos«  =  sin+sin40  +  cos4cos40cos ( A~A0 )  (2-60) 


Jhe  maximum  of  g  =  g(o)  with  respect  to  o  will  be  attained  at 
thus 


1& 

do 


«  *  If  ■  !£  -  = 

--  -  mss  +  3t,miKl  => 

0  *  *  3t*[0l  +  (2-61) 

which  vanishes  identically  for  o  s  0,  i.e.,  40  =  ♦,  A0  =  X.  This  is  the 
argument  originally  given  by  Bjerhammar  [1986,  p.  30].  However,  in 
order  for  o=0  to  be  a  maximum  (J*g/dwa)(o=0)  <  0  must  be  satisfied.  It 
holds  that 

0  *  h  t3t’sin"  ♦  *i~]  ■  3t*co.w+3t,(t*-l)  i;  (**§*)=> 

0  =  3t*coso  ♦  (2-62) 

Now  for  o  =  0  =>  si  no  =  0,  coso  =  1,  hence 
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d  =  [1  +  t*  -  2tcos«>]*  =>  d  =  1  -  t,  thus 
aatf  3t*(t*-ll 

^j5  («  =  0)  =  3t*  +  “YJ-tY*"'  or  l^>on 

0  <«  =  0)  =  (jilyx  tU-t)4  -  (l+t)]  (2-63) 

Now  for  0<t<l  =>  -l<-t<0  =>  0<l-t<l  =>  (l-t)4<l  (2-64) 

Also  0<t  =>  l<l+t,  thus 

(l-t)4<l<l+t  =>  (1-t)4  -  (l+t)<0  (2-66) 

which  yields  («=0)<0  and  thus  u=0  is  indeed  a  naxiaum  for  g. 

Therefore,  the  optimal  positions  for  the  Dirac  Impulses 
corresponding  to  gravity  anomaly  observations  are  at  the  nadir  points 
of  the  observation  stations. 

2.7  On  the  Optimal  Radius  of  the  Geosphere 

A  successful  application  of  the  Dirac  Impulse  method  requires  a 
suitable  choice  for  the  radius  of  the  internal  sphere.  Some 
suggestions  pertaining  to  gridded  data  can  be  found  in  KatBambalos 
[1981,  p.76]  and  Bjerhammar  [1986,  p.7]. 

A  suggestion  for  sparse  data  was  made  by  Sjoberg  [1978,  p.64] 
according  to  which 

r0  =  HE  -  RE  |  (2-66) 

where  RE  is  a  mean  earth  radius  (e.g.  RE=6371  km)  and  y  is  the 
average  angular  distance  between  neighboring  points.  Alternatively, 
one  can  attempt  to  compute  r0  from  the  given  data  if  this  is  possible. 
Two  methods  can  be  considered.  The  first  one  was  somewhat 
differently  regarded  in  Bjerhammar  [1986,  p.20).  He  considered  a 
least-squares  solution  whereas  here  a  minimum  norm  solution  will  be 
investigated.  The  system  in  (2-43)  can  be  written  as 

tm  =  a-A g*  =  F(X.)  (2-67) 

"ilh  *•  =  $*]•  *• 1  [if]-  *  *  [“f]  “d  v- 

Equation  (2-67)  represents  a  non-linear  system  in  the  unknowns. 
Linearising  by  Taylor's  theorem  and  neglecting  terms  of  order  2  and 
higher  yields 
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Aa  =  F(X<>)  +  (X.-Xo),  or  upon  letting 

»0  =  F(Xo);  t  =  ta-i0;  =  [^|f;  j*“]  =  B,  one  then  gets 

BX  -  i  =  0.  (2-68) 

In  the  case  of  n  observations  (2-68)  represents  a  system  of  n 
equations  in  (n+1)  unknowns  with  B  having  full  row  rank,  therefore 
the  condition  XTX=min  will  yield  (see  eq.  (2-49)  through  (2-55)) 

X  =  BT(BBT)*‘i  (2-69) 

with 
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2X  =  BT(BBT)-1Ej(BBT)-‘B.  (2-70) 

Now  let  us  write  B  as 

B  =  [G  a],  (2-71) 

where  a  =  r —  is  an  (nxl)  vector.  Denoting  by  ta  the  i-th  element  of 

■ 

d  la* 

ia,  the  i-th  element  a^  of  the  vector  a  will  be  aj  =  In  order 

to  evaluate  at  for  Ag,  (  and  q  observations  one  will  need  the  following 
derivatives 


...  at"  at  t" 

(l)  _  =  , 


(2-72) 


(ii)  iH>  =  is  iL  .  Recalling  (2-21)  and  (2-22)  one  gets 

3r0  3t  3r0  3r0  at  3r0 


t-cosw 


ad  _  t(t-coso) .  3u _ t_  ft 

3r0  “  r0d  ’  3r0  ~  2r0  l 

Now,  rewriting  the  i-th  equation  of  (2-67)  as 

4»,  =jll 

one  has 

*  I  U**  »«i. 

•«0  J=J  "”o  J 

where  for  gravity  anomalies,  from  (2-28) 


-  COScj 


')• 


(2-73) 


(2-74) 


(2-75) 


i 
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gif  =  ta^a)  -  3tscos«  -  ta 
and  for  vertical  deflectiona,  frc*  (2-34) 

$j)  •  *  (- -  h - 3-W ♦  —)— (S3- 


(2-76) 


(2-77) 


The  differentiation  of  the  g<  j  coefficients  required  in  (2-75)  yields 
(Appendix  A.3): 


=  £“[ss(1“2t*)  '  3-t^1~ta^t~COS°^  -  9tcos«  -  2] 


(2-78) 


for  gravity  anomalies.  For  vertical  deflections  it  yields  (Appendix  A.4) 


-b-3-^  ♦*—  <W- 


Equations  (2-78)  and  (2-79)  can  be  used  to  form  the  elements  of 
the  vector  a  in  (2-71).  Equation  (2-69)  can  be  used  to  solve  for  the 
vector  X  of  the  unknowns,  the  last  element  of  which  is  r0. 

For  the  actual  implementation  of  (2-69)  it  is  worth  noting  that 
since  G  in  (2-67)  is  of  full  rank,  then  [Boullion  and  Odell,  1971,  p.  18] 

GTl  -  CTlab 
b 


where  b  =  [1  +  (Cr*a)Ta-*a]-I(GT*a)TG-* 

Another  method  of  solving  for  r0  is  to  separate  the  data  in  two 
groups  and  of  n,  and  na  observations  respectively.  Denoting  by 
Ag*  the  Dirac  Impulses  corresponding  to  ij,  one  has  the  following 
linear  system: 

*  BAg*.  (2-80) 

Solving  (2-80)  for  Ag*  one  gets 

Ag*  = 

On  the  other  hand,  knowing  Ag*  one  can  predict  JJ  as 


(2-81) 
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*1  =  SAg* 

Substituting  (2-81)  in  (2-82)  one  gets 
A*  -  SR"li>  =  0 

Equation  (2-83)  can  serve  as  the  mathematical  model  of 
F(X„L,)  =  0  in  an  adjustment  scheme  with  r0  as  the  only 
Following  Uotila  [1985]  one  can  write 

X,  =  ,[r8]r,  X0  =  i[rg),;  X  =  1[6r0]1;  4,  =  [J-] 

(Ij+Oj  1 

^  =  [jg];  v  =  n  =  n‘+n* 

Linearizing  (2-83)  by  Taylor’s  theorem  and  neglecting  powers  of 
(X,-X0)  of  orders  2  and  higher  one  gets 

AX  +  BV  +  W  =  0  (2-85) 


(2-82) 


(2-83) 

the  form 
unknown. 


(2-84) 


with 

®  =  Ifs  =  I). 


(2-86) 


naxn 


naxnt  n2xn* 


Also 

A*  i?  ■  aS  ^  <*s ^  <*:>  -  W, 

naxl 

From  Dermanis  [1985,  p.187]  one  has  that 


d  /.«.  _  dA  B  .  dB  , 
dt  ~  *+*•  ®  +  *  **■  ’ 


dt 


dt 


—  1  =  -A-1  —  A-1 

dt  1  ~  A  dt  A  . 


(2-87) 


(2-88) 

(2-89) 


where  A,B  are  matrix-functions  of  the  variable  t  such  that  AB  and  A-1 
are  defined.  Applying  (2-88)  and  (2-89)  in  (2-87)  one  gets 


*  -  K1  &  '  Ifc)"-'! 


(2-90) 
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The  elements  of  the  matrices  A  and  B  are  given  by  (2-90)  and 
(2-86)  respectively.  The  least-squares  solution  of  (2-85)  is  given  by 
Uotila  [1985] 


P  =  *oaE?  (2-91) 

*b 

w  =  F(Xo,ih)  (2-92) 

M  =  BP-‘BT  (2-93) 

X  =  ~(ATM“lA)-lATM"lW  (2-94) 

X,  =  X0  +  X  (2-95) 

=  •.•(AH^A)"*.  (2-96) 


Since  the  problem  is  non-linear,  iterations  will  be  needed. 
Following  the  iteration  scheme  described  in  Uotila  [1982a]  one  can  write 
(for  the  i-th  iteration)  the  following: 


-  Evaluate  A,B  at  XJ  =  X*"1 ,  *A  =  i*m~l 

-  W»  =  F(XA,  <A)  +  B‘(4„-iA) 

-  m<  =  (b‘)p~1(b*)t;  p  *  9%  hi 

b 

X‘  =  -[(A<)t(M«)-‘A‘]_i(A*)7(M,)-1W< 
.  V,  =  -P~MB,)t(M*)-*(A,X‘  +  w1) 

|l«  =  l„  +  V* 

ixi  =  xA  +  x* 


(2-97) 

(2-98) 

(2-99) 

(2-100) 

(2-101) 


In  this  particular  case  it  is  worth  noting  that 

W  =  Ag  -  S}(R,)'in* 


(2-102) 


Yet  another  method  to  be  considered  for  optimal  r0  computations  is 
the  following.  A  measure  s*  of  the  quality  of  the  results  is  assumed  to 
be  a  second  order  polynomial  in  r0.  According  to  Bjerhammar,  [1986, 
p.17]  8a  is  the  variance  of  unit  weight  or  the  RMS  error  of  a  residual 
field  [ibid,  p.18]  (note  that  in  the  case  where  sa  =  RMS  (dig)  then  the 
units  of  b1  are  mgals  rather  than  mgal2).  Thus 


a*  =  a  +  br0  +  erg 


(2-103) 


The  value  of  r0  at  which  sa  is  an  extremua  will  be  the  root  of 
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j(i 

-r~  =  0,  or  equivalently  of  2cr0  +  b  =  0 
•ro 

The  root  of  this  equation  is 

r0  =  -  |^  (2-104) 

Furthermore,  for  c>0,  the  extremum  in  (2-104)  is  a  minimum. 
Assuming  that  three  different  radii  ro,,  roa,  ro3,  resulted  in  a},  si,  si 
respectively,  one  can  write  from  (2-103) 


s{  =  a  +  br0l  +  cr0 
s|  =  a  +  br0a  +  cr0 
sj  =  a  +  br0s  +  cr0 


(2-105) 


This  linear  system  of  three  equations  in  three  unknowns  a,b,c  can 
be  solved  to  yield  the  following 


a  =  s2 _ aj£glisLz.«JL  +  .r°M.r°.ii*i  ~  *i). 

(r0l-r0a)(r02-ro3)  (r0l-roa)(ro1-r0j 


b  _  (r«a+r«a)(sf-S?) 

(r,A+r0a)(,j-,i) 

(r*l-r°a)(r«J-r°l) 

(r°1~r«a)(r  oa~ro3) 

c  ■*  7  ■* 

fW 

e 

• 

i 

(r»1-r°J)(r«J-r°J) 

(fo.-ro^Cr.,-!-®,) 

(2-106) 


The  substitution  of  b  and  c  froa  (2-106)  in  (2-104)  yields 

.  (r»;-ro^(«i-«j)  -  (r.;-r,;)(«?-sj) 

0  ‘  2 ( (r9l-r0, ) (, j-, j )-(r0j-r0j) (a?_,g ) j  (Z"107) 

If  r«»  -  r0-h,  r«j  =  r0,  ro3  =  r«+h  is  selected,  one  obtains  from  (2-107) 


*•  ■  -  5fW-ST  <2-108> 

which  is  the  original  formula  derived  by  Bjerhammar  [1986,  p.18]  after 
correcting  a  minus  sign  error. 

The  actual  implementation  of  (2-107)  will  be  to  apply  the  predictors 
three  times  with  three  different  radii  and  record  the  resulting  s* 
values.  The  selection  of  the  three  initial  radii  is  quite  arbitrary. 
However,  in  order  to  make  (2-107)  most  effective  one  should  keep  in 
mind  that  (2-103)  is  a  parabola  and  its  minimum  will  be  best  computed 
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if  one  haa  a  point  in  the  ascending  part  of  the  curve,  one  point  in  the 
descending  part  of  the  curve  and  one  point  between  the  two. 


CHAPTER  III 


THE  BIHARMONIC  POTENTIAL  -  HARDY’S  METHOD 


3.1  Introduction 


The  gravitational  potential  of  a  mass  distribution  with  density  a 
occupying  a  volume  U  is  given  by  [Heiakanen  and  Moritz,  1967,  p.  3] 

V  =  GJ7/  %  do  (3-1} 

u  * 

where  G  =  Newton’s  gravitational  constant 
du=  differential  volume  element 
4  =  distance  between  du  and  the  evaluation  point. 

Let  the  volume  and  the  gravitational  potential  of  the  Earth  with 
mass  distribution  <rt  be  denoted  by  U ,  and  V,  respectively.  The 
gravity  potential  of  the  Earth  is  then  given  by 

W»  =  V,  +  ♦,  (3-2) 

where  ♦,  =  X«f(xa+ya)  is  the  centrifugal  potential  [ibid,  p.47)  and  o,  is 
the  Earth’s  angular  velocity. 

Introducing  the  standard  notion  of  the  reference  ellipsoid  with 
density  <r2,  volume  Ua,  gravitational  potential  Va,  gravity  potential  Wa 
and  angular  velocity  u  2  -o , ,  one  has 


v»  =  am  fj*  du 
u»  1 

(3-3) 

V,  =  am  f*  du 

«a  ‘ 

(3-4) 

Wa  =  v,  + 

(3-5) 

The  disturbing  potential  T  is  -fined  as 

T  =  W,  -  Wa. 

(3-6) 

Therefore 
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T  =W,-Wa  =V1+41-VJ-*J  =G///t1  do4  »i(x‘+y*)-aW%*  du-J-u^x’+y’) 

u.»  *  U-*  * 


and  since  &>t  =  «*a  one  gets 


T  =  am  du  -  G///  dM  +  G  /// 
«a  *  u*  u»_ua 


It  will  be  assumed  that  the  integral  of  over  the  difference  of 
the  two  volumes  Uj  and  Ua  is  negligible.  Denoting  by  <r  the  density 
«wn«uiiy  function,  i.e.,  a  -  9\ -aa  and  also  denoting  Ua  by  0  one  finally 

gets 


T  = 


dv 


(3-7) 


3.2  Hardy's  proposal  and  its  consaquencea 
3.2.1  Existence  of  the  biharmonic  potential 

The  representation  of  the  disturbing  potential  in  (3-7)  iB  singular 
at  points  that  induce  potential  since  at  these  points  1:0.  On  the  other 
hand,  since  there  are  infinitely  many  mass  distributions  that  have 
V,  as  their  exterior  potential  [Heiskanen  and  Moritz,  1967,  p.17],  there 
are  infinitely  many  density  anomaly  functions  <r  that  have  T  as  their 
exterior  disturbing  potential. 

Hardy  and  Nelson  [1986]  proposed  to  select  a  particular  family  of  <r 
functions,  namely  the  ones  that  are  zero  together  with  their  normal 
derivatives  at  the  boundary.  They  also  defined  a  function  p  such 
that 


p  =  lUv 


(3-8) 


where  the  Laplacian  operator  A  in  (3-8)  is  defined  as  usual,  i.e. , 


lz* 


and  it  refers  to  the  point  that  induces  potential.  Then  they  showed 
that  the  disturbing  potential  T  can  be  written  as 

T  =  G/j/pldv  (3-9) 

In  order  to  prove  (3-9)  one  needs  to  show  the  following:  Let  Q  be 
some  volume  of  mass  with  density  anomaly  function  a,  satisfying  a  s  0 
and  If/insO  at  the  boundary  S  of  0.  The  function  p  defined  by  2p=A<r 
is  such  that 
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X{J(7  '  pl)dv  =  0  (3_10) 

The  proof  will  in  principle  follow  Hardy  and  NelBon  [1986,  p.19].  One 
can  easily  establish  that  d(4/2)  =  1/4  at  points  with  4*0. 

We  will  distinguish  two  cases,  one  for  exterior  and  one  for  interior 
points,  due  to  the  singularity  at  4=0. 

la)  Bxterior  points: 

If  the  point  at  which  the  integral  in  (3-10)  iB  evaluated  is  outside 
the  volume  Q  then  by  Green's  identity  for  the  functions  4/2  and  a  one 
can  write  [Heiskanen  and  Moritz,  1967,  p.ll] 


A.  -  «A(f)j<!v  =  {/(73J 


• km*  -  {/is  k  -  f  $*■ 


gut,  on  S,  both  a  and  —  are  zero.  Therefore 

4n 

On  the  other  hand  p  =  and  d[^|  =  j,  hence 


jjf  iis.ailU  = 

s  12  3n  2  3nJ 


J^/(^  d<r  -  <rd(|]]dv  =  0  <=>  ~  j]dv  =  0,  q.e.d. 


(b)  Interior  points: 


Let  us  denote  by  0r  the  volume  that  will  remain  if  one  excises  from 
0  the  volume  of  a  sphere  with  center  at  P  and  radius  r.  The 
boundary  of  0r  is  denoted  by  Sr.  A  schematic  representation  is  given 
in  Figure  3. 


Figure  3. 


Derivation  of  the  Alternate  Integral  for  the  Disturbing 
Potential  at  points  Interior  to  the  Mass. 
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The  surface  ABCDA  is  the  boundary  S  of  0.  The  surface  ABCEA  is 
the  boundary  Sr  of  0r.  The  surface  ABC  is  SnSr  =  {points  QcS:  QcSr). 
The  surface  AEC  is  Sr-S  =  {points  Q«SP,  Q/S). 

By  Green's  identity  for  a  and  4/2  one  has 

v( i-  -  -d)]-  -  Fd is  - .  m*  -  "(§ is-fSS*- 

Mr  r  f* 

-."IfS-fUD*  (IS-iM)* 

The  first  integrr  l  in  the  above  formula  is  zero  since  a  -  3<r/3  n  =  0 
on  S  and  SnSrcS.  If  one  considers  Qr  small  enough  such  that 
(4/2)(4<r/4n)-(«'/2)(4J/an)  5  constant  =  c  [Scheik,  1986]  the  second 
integral  can  be  written  as  follows: 

(I  IS  -  f  4>  -  <=."  * *  c  W' 


which  of  course  goes  to  zero  like  ra.  Therefore 


tp (*  -  f)dv  *  -  *4(i))dv  1  -  “(I))*’] 

*  !2l"li  -  f  feW  *  °> 


r-*«lSr-S 


Therefore 


///(pi  -  f)dv  =  0  <=>  ///pidv  =  /£/  j  dv  <=>  G/^/pidv  =  Gf£f  |  dv 

(3-11) 

and  the  representations  of  T  in  (3-7)  and  (3-9)  are  identical. 


The  Poisson  equation 
AT  =  -4  ttG<r 


(3-12) 


corresponds  to  the  representation  (3-7)  for  the  disturbing  potential. 
The  same  way  as  above,  a  potential  represented  by  (3-9)  is  called 
biharmonic  because  then  it  satisfies 

+  &  +  B  +  +  2wh +  2^h - (3’13) 
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which  is  the  biharaonic  equation.  However,  the  derivatives 

a*9  a*9 

Jx1’  JyJ>  8za 


(3-14) 


must  exist  [Hardy  and  Nelson,  1986,  p.19].  Outside  the  masses,  where 
psO,  T  satisfies 

AaT  =  0.  (3-15) 

3.2.3  The  biharmonic  potential 

It  is  shown  in  Appendix  A.5  that  the  solutions  of  the  homogeneous 
biharmonic  equation  (AaT  =  0)  are 


T^r.O.X)  =•  I  r"  2  [(an-+racB-)cosaX+(bniB+r*dnm)sinaX]Pnw(cose) 

n-o  m =o 

(3-16) 

T.(r,0,X)=  l  -iqn-  2  [(an,+racn„)cos^+(bn«1+radnm)siruBX]Pnw(cos0) 

n=o  *  *=o 

(3-17) 

where  ann,  bnai,  cnM  and  dnM  are  arbitrary  constants.  Therefore  a 
biharmonic  function  T  can  be  represented  as 

T  =  Hi  +  raHa  (3-18) 

where  H,  and  Ha  harmonic.  If  Ha  vanishes  identically  then  T 
degenerates  to  a  harmonic  function. 


Equations  (3-16)  and  (3-17)  are  general.  Every  function  which  is 
biharmonic  inside  a  certain  sphere  can  be  expanded  into  a  series  (3-16) 
whereas  every  function  which  is  biharmonic  outside  a  certain  sphere 
can  be  expanded  into  a  series  (3-17). 

3.2.4  Further  consequences 

The  definition  pxXAe  is  a  partial  differential  equation  of  the  Poisson 
type.  Therefore,  the  above  definition  together  with  the  boundary 
conditions  a  :  89/8x1  -  0  uniquely  determine  a  [Kellogg,  1929,  p.215]. 
Actually,  since  the  earth  is  not  homogeneous,  one  should  write 


where  f,,  fa,  f3  are  functions  characterizing  the  inhomogeneity  of  the 
medium  [Volyskii  and  Bukhman,  1965,  p.36]. 

In  order  to  solve  p  s  XA 9  one  can  use  Green’s  third  identity 
[Heiakanen  and  Moritz,  1967,  pp.11-12]  to  get 
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JfJ  i  ».dv  =  -k.  ♦  //(i  £  -  .  |j(i))d. 


where 

k 


'4 n  if  the  evaluation  point  Q  is  inside  S, 

2 n  if  the  evaluation  point  Q  is  on  S, 

0  if  the  evaluation  point  Q  is  outside  S, 


and  n  is  the  nornal  to  S  directed  outward.  But  since  <r  =  —  •=  0  on 

*  an 

and  p  -  %A a  one  gets 


V 


£  dv  =  -  ^ 

a  '-*v  o 


(3-19) 


For  points  Q  inside  S  one  gets 


* 1  -  k  V  !  *•  (3-20) 

From  equation  (3-20)  one  can  see  that  the  singularity  of  T  at  4=0  is 
not  avoided.  It  is  simply  tranafered  to  a  singularity  in  a  at  the  same 
point  (4=0). 

An  obscure  point  in  Hardy’s  derivation  remains  the  existence  of  the 
fourth  order  partial  derivatives  of  T  in  (3-13).  The  reason  for  this  is 
that  at  least  one  of  the  second  order  partial  derivatives  of  T  must  be 
discontinuous  in  the  region  from  the  geocenter  to  infinity  [Heiskanen 
and  Moritz,  1967,  p.  S]  following  the  discontinuities  of  a. 

On  the  other  hand,  the  density  anomaly  function  9  is  assumed  to 
possess  properties  that  may  not  be  physically  reasonable.  At  first,  the 
second  partial  derivatives  of  a  are  assumed  to  exist.  Since  the  earth’s 
density  function  ia  very  likely  to  exhibit  discontinuities,  the  density 
function  of  the  reference  ellipsoid  must  be  discontinuous  in  such  a 
manner  that  both  9  and  its  partials  of  first  order  be  continuous. 
Furthermore,  the  density  function  of  the  reference  ellipsoid  muBt  be 
such,  that  9  together  with  its  normal  derivative  vanish  at  the 
boundary. 


The  aforementioned  requirements  of  the  method  are  not  justified 
from  the  point  of  view  of  the  physics  of  the  problem.  For  example, 
since  all  the  points  of  discontinuity  of  the  Earth’s  density  are  not 
known,  one  cannot  construct  a  reference  ellipsoid  such  that  the 
resulting  density  anomaly  function  9  and  its  partials  of  first  order  be 
continuous.  The  point  of  this  discussion  is  that  if  the  method  yields 
not  good  predictions  of  gravity  field  related  quantities,  this  should  come 
as  no  surprise  due  to  the  aforementioned  shortcomings  of  the  method. 
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3.3  Approximation  of  the  disturbing  potential 

Let  data  be  given  at  n  diacrete  points.  Let  us  also  subdivide  Q 
into  regions  Vt,  i=l,  2,  ...,  n.  The  biharmonic  sources  (sources  of 
biharmonic  potential)  are  defined  aB  follows  [Hardy  and  Nelson,  1986,  p. 
19] 


a,  =  ///  pdv 


(3-21) 


Let  P  be  the  point  at  which  the  disturbing  potential  will  be 
evaluated,  q(  be  one  of  the  given  data  points  and  q  be  any  point  in 
Vf.  Also,  let  *pq j  be  the  distance  from  p  to  q{,  4,  be  the  distance  of 
q  to  q,  and  4  be  the  distance  from  p  to  q. 
p,q,q{  can  be  written  as 


The  triangle  inequality  for 


1 4-4 


pqj 


d  4 


(3-22) 


Multiplying  both  sides  by  tpldv  and  fonaing  the  integral  for  V,  one 
gets 


(3-23) 


Summing  over  all  of  the  V,'s  one  gets 

gI  JJS  14-4  llpldv  *  G  ?  ///  4 , I p I dv  (3-24) 

t  =  l  V,  '  (  =  1 

Now 

J  2  G///4pdv  -  2  G///4pq  pdv j  =  g| JJJ4pdv  +. . .+/JJ4pdv-///4pqpdv 

»=i  v,  i=i  v,  v,  vn  vt 

-  '{"p^I  = 


*  G | III 4pdv  -  ///4pqipdv+. . .+  ///4pdv  -  ///4ptjnpdv j  < 

Vl  V1  Vn  Vn  " 

*  G|///4pdv  -  f/f4pSidvJ  +...+  g{ ///4pdv  -  i7J4p9npdv|  < 

VJ  vi  vn  Vn  " 

*  G///l4p-4piJi  J  •  lpldv+. . .+  G///j  4p-4piJn  j  •  Ipldv,  hence, 

Vl  Vn 


31 


I  |  G/JJipdv  -  I  G///ipq  pdvl  <  I  G///I lp-ipqi I- Ipldv 
1 1=1  v1  1=1  V,  1  1=1  Vf*  p  FH*' 

which,  upon  substitution  in  (3-24)  yields 

|g  2  fffi pdv  -  G  2  ///ipq  odvl  *  G  2  ///l,  Ipldv  (3-25) 

1  1  =  1  V,  1  =  1  V,  '  1  =  1  V, 

n 

Now  G  I  fffi pdv  =  T„  from  the  basic  integral  (3-9).  Also,  since  1D_ 

1=1  v,  F 

is  constant  for  each  V,  (3-25)  can  be  rewritten  as 

|Tp  -  G  2  ///pdv I  *  G  2  IfS*i  Ip'dv  (3-26) 

Let  us  denote  by  c«  the  distance  of  q«  to  the  furthermost  q  in  V(,  i.e., 
c |  >  *qq,  for  all  q  i  V|.  Also,  let  us  denote  by  e  the  maximum  e{,  i.e., 
t  -  max(tf)  for  i=l,2,...,n.  Then  i(  *  t  for  i=l,2,...,n  and  recalling 
(3-21),  (3-26)  can  be  written  as 

|tp  -  G  2  <pq  a,|  <  eG  2  ///  Ipldv  (3-27) 

Equation  (3-27)  implies  that  the  approximation 

T  =  G  I  i,a,  (3-28) 

i  =  i 

can  be  made  arbitrarily  good  by  an  appropriate  choice  of  t,  i.e.,  of  the 
size  of  the  subregions  V  (  (note  that  ipqi  was  substituted  by  I,  in 
(3-28)). 


3.4  Linear  functional  of  the  disturbing  potential 
3.4,1  Gravity  anomalies 

From  the  fundamental  equation  of  physical  geodesy  (2-7)  one  has 


Ag 


ar  _  2T 
Or  r 


(3-29) 


If  one  denotes  by  r  the  geocentric  distance  to  the  evaluation  point,  by 
r i  the  geocentric  distance  to  the  biharmonic  source  a(  and  by  +,X,4(,Xi 
their  geodetic  coordinates  respectively  one  has  from  (2-11)  and  (2-10) 

=  {ra+rf-2rrtcosw}11 


(3-30) 
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cos«  =  sin+sin+i  +  cos+cos+f cos(X-Xf ) . 


(3-31) 


I?  =  G,l1  “‘I?1'  With  H1  =  2iT  (2r-2r,co»«). 
4^  =  G  2  at  r~rlcos<J  ,  and  therefore 

*r  )  =  i 


Hence 


(3-32) 


3.4.2  Deflections  of  the  vertical 
Recall  equation  (2-30) 


ff]  =  i_  il  [coadl 

IqJ  yr  9u  Isinori  ' 


(3-33) 


The  required  derivative  is 


=  G  f  a<  =  G  2  at^rI'tain“,  thus,  using  (2-33)  one  gets 

dtt  )  =  i  ««  i  =  1  *t 

m  -  2  ?  a<r«  fcos+sin+j  -  sin4cos4<cos(X,-X)l  , 

IqJ  “  y  1=!  it  lcos4tsin(X,-X)  J*  ' 


(3-34) 


If  one  places  the  biharmonic  sources  on  the  surface  of  the 
geosphere  with  radius  r0,  equations  (3-32),  (3-34)  will  become 


*■-■!.  ■1[E¥2SS*i£] 


ffl  _  G  8  a<rn  fcostsin+t  -  si 
W  ~  T«=i  1  lcos41sin(X<-X) 


-  sin+cos+fCos(Xt~X] 


(3-35) 


(3-36) 


I  *  {ra+rj  -  2rr0cos«}* 
Introducing  t  and  d  by 


(3-37) 


i 
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and  d  =  {l+ta-2tcos«}* 
one  gets 

1  =  {ra+rj-2rr0cos«»},t  =  r{l+ta-2tcosu}*  =  rd;  hence 

U  -  -8,1, 

/(]  _  §t  B  fcos+sin4(  -  sin4cos4tCos(X,-X) 1 
y  d  lcoa4tsin(Xi-*)  > ' 

Finally,  letting 

c,  =  Ga, 

one  can  write  (3-40),  (3-41)  as 


(3-38) 

(3-39) 


(3-40) 

(3-41) 


(3-42) 


Ag  =  +  2d]c<$  (3-43) 

m  _  t  a  cx  fcostsin+i  -  sin4cos4,cos(X,-X)|  (3-44) 

li|i  y  t=t  d  lcos4|Sin(X(-X)  >  ‘ 

Equations  (3-43),  (3-44)  can  be  used  to  compute  c*  from  observed 
Ag,  (  and  7i  and/or  to  predict  Ag,  (,  v  from  previously  computed  c, 
values.  The  associated  linear  system  and  its  solution  will  be  identical 
to  the  one  described  in  Section  2.4.  Also,  the  propagation  of  the  data 
noise  into  the  predictions  will  be  performed  in  a  manner  identical  to 
the  one  described  in  Section  2.5. 


CHAPTER  IV 


THE  DATA 


4.1  Introduction 


The  White  Sands  Teat  Area  is  located  at  the  western  outskirts  of 
the  Rio  Grande  Rift  System  in  New  Mexico.  The  tests  of  the  two 
methods  were  performed  with  data  in  the  portion  of  White  Sands  bound 
by  the  parallels  32*N  and  34*N  and  the  meridians  253*E  and  254*E. 
This  area  is  mainly  a  plateau  at  a  level  of  1200  m  to  1400  m  (Figure  4). 
The  Oacura  and  San  Andres  mountain  chains  cross  the  area  in  a 
North-South  direction.  The  geological  constitution  of  the  area  is  mainly 
young  meaozoic  sediments  complimented  by  some  late  tertiary  volcanics 
[Schwarz,  1983,  p.  21. 

The  bulk  of  the  White  Sands  Test  data  were  made  available  to  the 
Special  Study  Group  4.70  of  the  International  Association  of  Geodesy  by 
the  National  Geodetic  Survey,  NOS,  NOAA,  Rockville,  Maryland.  C.C. 
Tscherning  did  some  initial  data  screening  and  then  arranged  the 
different  files  for  the  tests  [Schwarz,  1983].  The  test  data  used  were 
made  available  to  us  by  C.C.  Tscherning  and  were  identical  to  the  data 
used  in  the  collocation  solution  of  the  report  "White  Sands  Revisited" 
[Kearsley  et  al.,  1985]. 


4.2  The  Two  Solutions 

For  each  method  two  independent  solutions  were  employed.  One 
for  the  (l#xl*)  area  bound  by  parallels  33*N  and  34*N  and  called  the 
North  Block  (NB)  and  one  for  the  (l*xl*)  area  bound  by  parallels  32*N 
and  33*N  and  called  the  South  Block  (SB). 
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Figure  4.  Topographic  Map  of  the  New  Mexico  Test  Area  From  a  2  x2' 
DTM  (Cl  =  100m). 
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4.3  Gravity  Anomalies 

The  gravity  anomaly  data  are  free-air  values  referenced  to  the 
Geodetic  Reference  System  1967.  Their  geodetic  latitude  and  longitude 
are  given  in  the  NAD27,  their  height  is  orthometric  and  their  standard 
deviation  is  *  2  mgals.  There  are  384  observations  in  the  NB  (Figure 
5),  548  observations  in  the  SB  (Figure  6),  82  control  values  in  the  NB 
(Figure  7)  and  123  control  values  in  the  SB  (Figure  8). 


4.4  Deflections  of  the  Vertical 

The  vertical  deflection  data  are  astrogeodetic  values  referenced  to 
the  NAD27.  Their  geodetic  coordinates  are  given  in  the  NAD27,  their 
height  is  orthometric  and  their  standard  deviation  is  *0.'3.  There  are 
67  (£tv)  observed  pairs  in  the  NB  (Figure  9),  63  observed  pairs  in  the 
SB  (Figure  10),  176  control  pairs  in  the  NB  (Figure  11)  and  208  control 
pairs  in  the  SB  (Figure  12). 

4.5  Conversion  of  the  data  to  an  approximately  geocentric  system 

The  system  in  which  all  the  calculations  were  carried  out  was  an 
approximately  geocentric  system  with  the  ellipsoidal  parameters  of 
GRS80.  The  datum  transformation  parameters  from  NAD27  to  the  new 
system  are  [Schwarz,  1983,  p.  13],  [Tscherning,  1987] 

Ax  =  -22m,  Ay  =  157m,  Az  s  176m,  *  =  0,  i>  -  0,  u  -  -017,  AL  =  0  (4-1) 


The  geodetic  latitude  and  longitude  can  be  transformed  to  the  new 
system  by 


♦NA0a7+d4j 

^NADa7+d^' 


with  [Rapp,  1981,  pp. 70,77] 


(4-2) 


sin+cosA  .  sin+sinA  .  .  cos4  .  easin4cos4 

"  ’  -  — SJh  *x  -  TCh *»  +  iwT  iz  *  W(M+h) 

sin4cos4(2N+e,aMBina4) 

Mfh  U 

sinA  .  .  cosA  , 

(N+h)cos4  X  (N+h)cos4  y  “ 


and  [Rapp,  1984,  pp. 21, 30, 35] 
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Figure  6.  Distribution  of  the  548  Gravity  Observations  at  the  South 
Block  of  the  White  Sands  Test  Area. 
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Figure  9. 
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Distribution  of  the  67  Observed  Vertical  Deflection  Pairs  at 
the  North  Block  of  the  White  Sands  Test  Area. 
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Figure  10.  Distribution  of  the  63  Observed  Vertical  Deflection  Pairs  at 
the  South  Block  of  the  White  Sands  Test  Area. 


W*  =  l-e*sina4;  M  = 


(4-4) 


Also  [Rapp,  1984,  p. 169] 

&NA037  ~  6378206.4  a 
fNAD27  =  1/294.978698 
^neh  =  6378137  ■ 
fMEH  =  1/298.257222101 

and  Aa  =  aNEH  -  aMAD17 

~  ^MEH  “  ^NA037 

The  error  in  using  orthometric  height  instead  of  geometric  height  in 
(4-3)  is  less  than  O'.OOl. 

Similarly  for  the  vertical  deflections  one  has 

[(men  =  f«AD27  +  d$1  f4-5^ 

Uses  -  (nao27  +  d^J 

where  [Rapp,  1981,  p.  74} 

(2 ;  :2LJ  <4-e> 

The  relation  between  normal  gravity  computed  with  the  GRS80  and 
the  GRS67  reference  ellipsoid  is  [Schwarz,  1983,  p.  13] 

7i ,go  =  7*,.7  +  (0.8316  +  0.0782sin24  -  0.0007sin*4)  (4-7) 

Furthermore,  in  Section  2.2.1  it  wbb  assumed  that  the  space  outside 
the  boundary  is  empty  which  implies  that  the  atmospheric  and  the  tidal 
effects  have  been  removed  from  the  observed  gravity.  As  far  as  the 
tidal  corrections  are  concerned  it  is  assumed  that  they  have  been 
modeled  during  the  observation  reduction  process.  The  atmospheric 
corrections  will  be  computed  by  [Wichiencharoen,  1982,  p.  5] 

«gA  =  (0.8658  -  9.727xlO-5H  +  3.482xlO-*H*)agals  (4-8) 

where  H  is  the  orthoaetric  height  in  aeters. 

Therefore  the  gravity  anomalies  referenced  to  GRS80  and  corrected 
for  atmospheric  effects  are  given  by 

Agcsiso  ~  Agg r $ « 7  —  0.8316— 0.0782sin24NEM  +  0.0007sin44nEn  +  ^gA 

(4-9) 

Contour  maps  of  the  observations  at  both  the  North  and  the  South 
Block  of  the  White  Sands  Test  Area  are  shown  in  Figures  13  through 
18. 
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Figure  13.  Contour  Map  From  the  Observed  Gravity  Anomalies  at  the 
North  Block  of  the  White  Sands  Test  Area.  (CI=20  mgals). 
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Figure  14.  Contour  Map  From  the  Observed  Meridional  Deflections  at 
the  North  Block  of  the  White  Sands  Test  Area.  (01=2”). 
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Figure  IS.  Contour  Map  From  the  Observed  Prime  Vertical  Deflections  at 
the  North  Block  of  the  White  Sands  Test  Area.  (CI=2"). 


Figure  16.  Contour  Map  From  the  Observed  Gravity  Anomalies  at  the 
South  Block  of  the  White  Sands  Test  Area.  (CI=20  mgals). 
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4.6  Removal  of  Reference  Field  and  Residual  Terrain  Model  (RTM) 
Effects 


The  predictors  will  be  applied  to  the  mid-frequencies  of  the 
anomalous  potential  spectrum.  The  low  frequency  effects, 

corresponding  to  wavelengths  of  about  1*  (~111  km)  or  larger  will  be 
accounted  for  by  the  OSU86F  reference  field  to  degree  and  order  180 
[Rapp  and  Cruz,  1986].  The  reference  field  computations  will  be 
performed  as  described  in  [Rapp,  1982b]. 

The  Residual  Terrain  Model  (RTM)  effect  is  the  effect  of  the  masses 
between  the  actual  topography  and  a  mean  elevation  surface,  on  the 
gravity  anomalies  and  the  deflections  of  the  vertical.  The  actual 
topography  is  represented  by  a  detailed  OTM  whereas  the 
aforementioned  mean  elevation  surface,  termed  the  reference  surface,  is 
a  coarse  DTM  which  is  usually  derived  from  the  detailed  DTM  by 
averaging.  The  residual  topography  is  modeled  as  an  assembly  of 
rectangular  prisms  with  a  constant  positive  or  negative  density 
depending  on  whether  the  terrain  surface  is  above  or  below  the 
reference  surface  [Kearsley  et  al.,  1985,  p.  53].  The  effect  of  the 
residual  topography  on  the  gravity  anomalies  and  the  vertical 
deflections  is  removed  computationally,  so  that  the  residual  quantities 
refer  to  the  reference  surface  rather  than  the  actual  topography 

[Forsberg,  1988].  Therefore,  the  RTM  effects  account  for  the  short 
wavelength  features  of  the  anomalous  potential  spectrum  [Kearsley  et 
al.,  1985,  p.  55]. 

The  question  of  the  optimal  RTM  computations  has  not  been 
completely  answered  yet.  For  example,  Forsberg  and  Tscherning  [1981] 
used  two  different  grid  sizeB  as  reference  surfaces  for  testing 
purposes.  On  the  other  hand,  Cruz  [1985,  p.  74]  used  a  spherical 
harmonic  expansion  of  the  topography  as  the  reference  surface. 

Moreover,  Kearsley  et  al.  [1985,  p.  55]  performed  tests  using  different 
DTMs  as  reference  surfaces.  These  tests  indicated  that  the  coarser  the 
reference  surface  the  smoother  the  residual  field.  However,  different 
reference  surfaces  have  an  insignificant  impact  on  the  predictions  due 
to  the  remove-restore  operation  [ibid,  p.  55]. 

The  RTM  effects  used  in  this  investigation  were  computed  by  R. 
Forsberg  and  C.C.  Tscherning  [Forsberg,  1988].  The  same  RTM  effect 
values  were  also  used  in  other  tests  with  the  White  Sands  data 

[Schwarz,  1983]  and  [Kearsley  et  al.,  1985].  A  (30"x30")  elevation  grid 
which  extends  over  the  area  31*30'  <  ♦  <  35*  r.nd  252*  <  X  <  255*  was 
used  as  the  detailed  DTM,  whereas  a  (30'x30')  grid  was  used  as  the 
reference  surface.  Consequently,  the  majority  of  the  signal  of  the 
anomalous  field  at  wavelengths  of  30'  {'*'55  km)  or  smaller  was  removed 
by  the  RTM  computations.  The  remaining  part  of  the  signal  (between 
wavelengths  55  km  and  110  km)  was  left  to  be  handled  by  the 
predictors  [Forsberg,  1988]. 
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Tables  1  and  2  show  the  results  of  the  removal  of  both  the 
reference  field  and  the  RTU  effects  from  the  observed  and  the  control 
values  respectively.  The  residual  quantities  V  in  these  tables  are 
defined  as 


VAg  =  Agcnsao  ~  4«osu»«f  ~  a«rtm] 

V$  =  (sen  ~  loSU««F  “  fllT*  (4-10) 

Vij  *  *?nem  ~  vo%vm»r  -  nsTs  * 


Table  1.  RMS  Values  of  OSU86F  and  RTM  Effects  on  Observations. 


OBSERVATION 

E3HE53 

Ag(agal) 

26.59 

14.86 

19.64 

17.68 

e  <") 

3.73 

U'M.'m 

1.95 

1.80 

..a  in _ 

6.42 

4.53 

3.34 

4.09 

Table  2.  RMS  Values  of  OSU86F  and  RTM  Effects  on  Control  Data. 


1 

RTM 

LiMiiliMUl: 

25.54 

10.94 

15.95 

19.35 

■lHM 

3.82 

2.66 

1.76 

1.63 

■m 

6.78 

4.65 

3.42 

4.33 

Tables  1  and  2  indicate  that  the  residual  field  is  smoother  than  the 
original  field.  Figures  19  through  24  show  contours  from  the  residual 
observations  in  both  the  NB  and  the  SB.  Comparison  of  Figures  13  to 
19,  14  to  20t  15  to  21  for  the  North  Block  and  16  to  22,  17  to  23,  18  to 
24  for  the  South  Block  reveals  that  the  basic  signature  of  the 
anomalous  potential  is  not  lost  by  the  removal  of  the  OSU86F  and  RTM 
effects.  However,  some  irregularities  of  the  original  field  have  been 
smoothed  out  by  these  computations. 

The  CPU  time  required  to  compute  reference  field  effects  is  about 
0.5  aec/station  and  to  compute  RTM  effects  is  about  0.5  sec/station  on 
the  IBM  3081  main  frame. 

In  conclusion,  the  computation  scheme  will  be  to  remove  the 
OSU86F  and  RTM  effects  from  the  observations,  perform  the  predictions 
with  the  residual  field  and  then  restore  the  OSU86F  and  RTM  effects  at 
the  prediction  stations. 
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Contour  Map  of  the  Residual  Observations  at  the  North 
Block  of  the  White  Sands  Test  Area.  Gravity  Anomalies 
(CI=10  ragals). 


Figure  20. 


Contour  Maps  of  the  Residual  Observations  at  the  North 
Block  of  the  White  Sands  Test  Area.  Meridional  Deflections 
(Cl=l"). 
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Figure  21.  Contour  Map  of  the  Residual  Observations  at  the  North 
Block  of  the  White  Sands  Test  Area.  Prime  Vertical 
Deflections  (CI=1">. 
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Figure  22.  Contour  Map  of  the  Residual  Observations  at  the  South 
Block  of  the  White  Sands  Test  Area.  Gravity  Anomalies 
(CI=10  mgals). 


CHAPTER  V 


ANALYSIS  OP  THE  RESULTS 


5.1  Introduction 


The  Bjerhammar  and  Hardy  predictors  were  tested  with  data  in  the 
portion  of  the  White  Sands  Test  Area  bound  by  parallels  32*N  and  34*N 
and  meridians  253 *E  and  254  *E.  For  each  method  tested  two 
independent  solutions  were  employed.  One  for  the  l*xl*  area  bound  by 
parallels  33*N  and  34*N  and  called  the  North  Block  (NB)  solution  and  one 
for  the  l*xl*  area  bound  by  parallels  32*N  and  33*N  and  called  the 
South  Block  (SB)  solution. 

Two  factors  that  contribute  greatly  to  the  quality  of  predictions  of 
quantities  related  to  the  Earth’s  gravity  field  are  the  topography  and 
the  data  coverage  of  the  area  of  interest.  The  2*xl*  area  at  which  the 
two  predictors  were  tested  presents  significant  variations  in  both  of 
these  factors.  As  far  as  the  terrain  is  concerned,  Figure  4  reveals 
relatively  large  valleys  and  rather  high  mountain  peaks,  the  NB  being 
more  irregular  than  the  SB.  In  relation  to  data  coverage,  the  SB  is 
superior  to  the  NB  in  terms  of  gravity  anomalies.  It  possesses  more 
observations  the  distribution  of  which  is  more  even  than  the  ones  of  the 
NB  (Figures  5  and  7).  However,  in  terms  of  vertical  deflections,  Figures 
6  and  8  show  that  the  NB  is  superior  to  the  SB  as  related  to  both 
amount  and  distribution  of  data. 

Furthermore,  both  the  terrain  and  the  data  coverage  vary 
significantly  within  blocks.  Therefore,  in  order  to  understand  the 
performance  of  each  method  better,  the  area  was  divided  into  eight 
0'.5x0'.5  sub-blocks  (hereafter  to  be  referred  to  as  "the  O’. 5  blocks”  or 
simply  "blocks”)  following  Kearsley  et  al.  [1985,  p.61].  The  eight  0:5 
blocks  are  shown  in  Figure  25. 
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Figure  25.  The  eight  blocks  used  to  test  the  performance  of  the 
two  predictors  at  the  White  Sands  Test  Area. 


The  "goodness  of  representation"  factor  R  can  be  defined  as  follows 
[Kearsley  et  al.,  1985,  p.63] 

R  =  10-3C*„  (5-1) 

where  <rH  is  the  RMS  height  variation  with  respect  to  the  mean  elevation 
of  the  area  and  C  in  kma/station  is  the  average  coverage  of  the  area 
defined  as 

C  =  jj  (5-2) 

where  A  is  the  area  and  n  is  the  number  of  stations  in  A. 

A  large  R  value  indicates  either  few  observations  or  highly  varying 
terrain  or  both.  It  is  essential  to  realize  that  R  is  a  relative  quantity 
for  intercompari8on  of  the  eight  0*.5  blocks.  This  is  to  say  that  an  R 
value  of  10  for  the  observed  deflections  represents  a  good  sub-block 
whereas  it  may  represent  a  very  poor  sub-block  in  terms  of  gravity 
control  stations.  As  a  specific  example,  block  #8  with  R&g  =  6.1  for  the 
control  stations  is  considered  as  representing  the  area  well  whereas 
blocks  *2  and  3  with  R^g  equal  to  5.7  and  5.3  respectively  for  the 
observation  stations  is  considered  as  representing  the  area  poorly.  The 
following  table  shows  these  details  by  sub-block.  In  this  table  P  stands 
for  Poor,  M  for  Medium  and  G  for  Good  sub-blocks. 


Table  3.  Terrain  Characteristics  and  Data  Coverage  in  the  Test 
Area. 


TERRAIN 

TYPE 


Hilly  +  Flat 
Mountainous 
Mountainous 
Hilly  ->  Flat 
Mountainous 
Hilly  ->  Flat 
Hilly  ->  Flat 
Hilly  ->  Flat 


OBSERVATIONS 


R  TYPE 


CONTROL  STATIONS 


81  25 
200  29 
217  24 
128  30 
176  19  173 
22  144 
187(18  238 
72  18  138 


M 
P 
P 

18.4  M 

30.5  M 
13.0  Q 
44.4  M 


69  114 
84  39 


M 

M 

P 

M 

M 

M 

M 

M 

M 

M 

M 

G 

M 

M 

G_J 

G 

Table  3  quantifies  the  differences  by  sub-block  in  terms  of  both 
terrain  characteristics  and  data  coverage.  For  example,  the  difference 
between  the  SE  30'x30'  portion  of  the  SB  (sub-block  #7)  and  the  SW 
30'  x30'  portion  of  the  SB  (sub-block  #8)  in  terms  of  vertical  deflection 
observations  is  clearly  demonstrated  by  the  R^>7)  factors  which  are  44.4 
and  9.9  respectively.  From  Table  3  one  can  see  that  the  SB  solution 
should  be  considered  more  representative  of  the  capabilities  of  the 
method  under  consideration  than  the  MB  one.  Furthermore,  in  terms  of 
gravity  anomaly  observations,  sub-blocks  #2  and  3  are  not  anticipated  to 
contribute  greatly  to  a  possible  good  NB  solution,  whereas  sub-blocks  #6 
and  8  are  capable  of  being  the  major  contributors  to  a  possible  high 
quality  of  the  SB  solution.  On  the  other  hand,  the  gravity  anomaly 
control  data  coverage  in  sub-block  *2  is  poor  rendering  the  comparison 
of  results  not  very  reliable  in  this  block,  whereas  block  #8  is  good  for 
this  purpose.  Also,  reliable  sub-blocks  for  vertical  deflection 
comparison  in  terms  of  control  values  are  blocks  6  and  8. 


In  order  to  evaluate  the  two  predictors  the  differences 


dAg  =  Agc  -  Ag. 

<4  =  (c  -  fp 

dn  =  Vc  ~  *?p 


(5-3) 


will  be  examined  where  the  subscript  c  refers  to  the  control  values  and 
the  subscript  p  refers  to  the  results  of  either  the  Bjerhammar  or  the 
Hardy  method.  The  differences  in  (5-3)  are  due  to  errors  in  the 
prediction  as  well  as  errors  in  the  control  data. 


One  of  the  moat  important  factors  influencing  the  quality  of  the 
predictions  with  the  Bjerhammar  method  is  the  radius  r0  of  the  internal 
sphere.  Up  to  a  certain  extent  r0  is  a  coupling  factor  in  the  sense  that 
the  improvement  of  the  predictions  expected  by  a  smooth  terrain,  by 
good  data  coverage  and  by  the  removal  of  reference  field  and  residual 
terrain  model  effects  can  be  easily  nullified  by  an  unsuccessful  choice 
of  r0.  More  importantly,  an  inappropriate  choice  of  rQ  can  render  the 
downward  continuation  impossible.  Due  to  the  aforementioned  effects  of 
r0  on  the  predictions,  efforts  were  made  to  compute  it  from  the  data. 

Two  methods  never  before  tested  with  the  Bjerhammar  predictor 
were  attempted.  The  first  method  is  the  minimum  norm  (pseudo)  solution 
given  by  equations  (2-67)  to  (2-70).  It  was  tested  in  the  NB  with  384 
gravity  anomalies  as  observations.  The  unknowns  were  both  384  Dirac 
Impulses  and  the  optimal  radius  of  the  geosphere,  a  total  of  385 
unknowns.  An  approximate  value  of  6350  km  for  the  optimal  radius  r0 
resulted  after  2  iterations  in  an  adjusted  value  of  6350  *  0.17xl0~s  km. 
The  residuals  were  in  the  order  of  10~9  mgals,  VTPV  was  7xl0-15  and 
the  standard  deviations  of  the  Dirac  Impulses  exceeded  the  values  of  the 
Impulses.  Similar  results  were  attained  after  two  iterations  when  the 
approximate  value  for  r0  was  6360  km.  The  order  of  magnitude  of  the 
residuals  and  VTPV  can  be  explained  by  the  fact  that  no  redundant 
observations  are  present  in  the  solution.  The  large  standard  deviations 
of  the  Impulses  stress  that  the  values  for  the  Impulses  are  evaluated 
with  very  large  uncertainty.  The  only  peculiar  result  is  the  small 
standard  deviation  of  r0  even  though  the  adjusted  value  of  r0  is  the 
same  one  as  the  approximate.  At  any  rate  this  method  did  not  seem  to 
have  computed  an  optimal  radius  of  the  geosphere. 

The  second  method  to  compute  r0  from  the  data  is  to  separate  the 
observations  into  two  groups  and  to  consider  the  first  group  as 
observed  values  and  the  second  one  as  control  values.  This  type  of 
solution  is  given  by  equations  (2-80)  through  (2-102).  In  this  case  we 
only  have  one  unknown,  namely  the  optimal  geosphere  radius  r0  and 
therefore  we  have  only  one  normal  equation.  This  method  was  tested  in 
the  NB.  The  results  of  the  solution  (iteration  #0)  as  well  as  some 
selected  iterations  are  shown  in  Table  4.  In  Table  4,  N  is  the  normal 
matrix  (of  dimensions  (lxl))  and  U  is  the  right  hand  side  vector  (of 
dimensions  (lxl))  of  the  normal  equations. 
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Table  4.  Data  Separation  Method  of  Computing  r0.  Bjerhammar 
Predictor. 


ITERATION  « 

E7W251 

N 

U 

Arn(km) 

rg  (km) 

0 

6360.000 

0.000  015  520 

0.045  439  743 

-2.928 

6357.072 

1 

6357.072 

0.000  007  390 

0.031  566  632 

-4.272 

6352.801 

5 

6346.544 

0.000  001  266 

0.002  313  902 

-1.827 

6344.717 

10 

6339.725 

0.000  000  799 

0.000  635  825 

-0.816 

6338.909 

20 

6335.527 

0.000  000  578 

0.000  082  876 

-0. 143 

6335.384 

30 

6334.783 

0.000  000  546 

0.000  019  254 

-0.035 

6334.748 

40 

6334.576 

[>Hi  •  iMl  • 

0.000  003  724 

-0.007 

6334.569 

45 

6334.549 

[till  i  iV  t 

0.000  001  777 

-0.003 

6334.546 

50 

6334.536 

0.000  000  535 

0.000  000  948 

-0.002 

6334.534 

52 

6334.533 

0.000  000  535 

0.000  000  767 

-0.001 

6334.531 

53 

6334.531 

0.000  000  535 

0.000  000  485 

-0.001 

6334.530 

From  Table  4  one  observea  that  the  iteration  criterion,  the 
correction  to  r0  be  leas  than  1  m,  was  met  after  53  iterations. 
Furthermore,  the  normal  matrix  stabilizes  only  at  the  45th  iteration  and 
the  rate  at  which  the  correction  to  r0  tends  to  zero  is  very  low.  Most 
importantly,  the  resulting  adjusted  value  of  rg  =  6334.530  km  yields  RMS 
differences  of  control  minus  predicted  values  in  the  order  of  7.73  mgals 
for  Ag,  2CT.'53  for  (  and  30107  for  rj.  These  differences  are  much  larger 
than  the  ones  yielded  with  the  same  data  type  when  an  optimal  radius 
was  computed  prior  to  the  solution  as  it  will  be  demonstrated  in 
Subsection  5.2.2.I. 

The  overall  conclusion  from  both  of  the  aforementioned  methods  is 
that  they  did  not  yield  an  optimal  geosphere  radius.  Therefore  the  s* 
method  given  by  (2-103)  to  (2-108)  will  be  used  for  optimal  geosphere 
radius  computations. 


5.2.2  The  Asymmetric  Kernel  Approach 

The  predictor  defined  by  equations  (2-8)  through  (2-34)  is  called 
the  Asymmetric  Kernel  (AK)  approach  to  be  distinguished  from  the 
Symmetric  Kernel  (SK)  Approach  given  by  (2-35)  through  (2-42).  The 
Symmetric  Kernel  Approach  is  given  its  nama  by  Bjerhammar  [1986,  p. 
48).  In  the  SK  approach,  t  in  (2-38)  is  a  symmetric  quantity,  i.e., 
invariant  with  respect  to  i  and  j.  Recall  that  t  of  the  AK  is  not 
symmetric  (  see  equation  (2-15)).  Furthermore,  if  the  observations  are 
only  gravity  anomalies,  then  the  design  matrix  G  in  (2-43)  is  symmetric. 
The  Asymmetric  Kernel  Approach  was  given  this  name  in  this 
investigation  in  order  to  distinguish  it  from  the  SK  Approach.  In  what 
follows  the  Dirac  Impulses  will  be  located  at  the  nadir  points  of  the 
observations.  In  the  tables  that  follow  the  differences  in  Ag,  (  and  tj 
are  control  minus  predicted.  The  notation  r0A9,  r0f  and  r0 will  be 
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differences  in  Ag,  (  and  tj  respectively  as  sJ  values.  On  the  other 
hand,  r0  will  be  defined  as: 

r0  =  5  (rcA9  +  r0f  +  To7))  (5-4) 

Every  variation  of  each  predictor  will  be  tested  and  compared  with 
similar  results  after  the  following  cycle  is  completed: 

(1)  Perform  three  solutions a  with  three  different  radii, 

(2)  Compute  r049,  r0$  and  r0V  by  (2-107), 

(3)  Compute  r0  by  (5-4). 

(4)  Perform  the  final  solution  with  r0=r0. 

5.2.2. 1  Prediction  Using  Only  Gravity  Anomaly  Data 

In  the  case  where  only  gravity  anomalies  are  observed  the  exact 
solution  as  given  by  (2-44)  and  (2-45)  applies.  The  elements  of  the 
design  matrix  G  in  (2-44)  and  (2-45)  are  given  by  (2-28).  The  results 
of  both  the  NB  and  the  SB  solutions  with  three  different  radii  are 
shown  in  Table  5. 


Table  5.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Asymmetric  Kernel  Approach  and  Only  Ag 
Observed.  Bjerhammar  Method. 


SOLUTION 

r0(kn) 

Ag(mgal) 

*n 

V  n 

6355 

4.03 

0.99 

0.99 

wortn 

6360 

3.32 

0.91 

0.97 

dIOCK 

6365 

2.96 

0.92 

1.11 

Caulk 

6355 

3.87 

1.00 

1.17 

oouin 

6360 

3.80 

0.92 

1.13 

d  iOCK 

6365 

3.89 

0.89 

1.19 

Using  the  results  of  Table  ^5  with  the  s*  method  (eq.  (2-107))  and 
the  RMS  differences  one  gets  r0  -  6362.571  km  for  the  NB  and  r0  s 
6361.562  km  for  the  SB  solution.  The  results  of  the  NB  and  the  SB 
solutions  with  r0  =  r0  are  shown  in  Table  6. 
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Table  6.  RMS  Differences  Between  Predicted  and  Control  Values  with 
the  Asymmetric  Kernel,  Only  Ag  Observed  and  the  Optimal 
Radius  of  the  Geosphere  r0,  by  015  Block.  Bjerhammar 
Method. 


I  BLOCK  1  BMS  DIFFERENCES  1 


North 

3.08 

0.90 

1.00 

1 

4.27 

0.91 

0.81 

2 

2.54 

0.64 

0.99 

3 

2.84 

1.26 

0.70 

4 

2.86 

0.81 

1.28 

South 

3.79 

0.90 

1.13 

5 

3.43 

1.27 

1.36 

6 

6.22 

0.94 

1.28 

7 

3.20 

0.58 

0.75 

8 

2.23 

0.74 

0.89 

From  Table  6  one  observes  that  the  best  Ag  predictions  were 
performed  at  sub-block  18  and  the  best  (  and  tj  predictions  were 
performed  at  sub-block  17. 

5.2.2.2  Prediction  Using  Both  Gravity  and  Vertical  Deflection  Data 

In  the  case  where  both  Ag  and  (f,q)  are  observed  the  least  squares 
solution  applies  with  as  many  degrees  of  freedom  as  deflection  pairs. 
This  type  of  solution  is  given  by  (2-46)  through  (2-48)  and  the  elements 
of  the  design  matrix  G  are  given  by  (2-28)  for  gravity  anomaly 
observations  and  by  (2-32)  for  vertical  deflection  observations.  The 
results  for  three  different  radii  are  given  in  Table  7. 


Table  7.  RMS  Differences  Between  Predicted  and  Control  Values  with 
the  Asymmetric  Kernel  Approach  and  Both  Ag  and  (f,»j) 
Observed.  Bjerhammar  Method. 


tn 

Laill 

North 

Block 

6355 

6360 

6365 

3.78 

3.56 

3.76 

0.87 

0.78 

0.77 

0.73 

0.73 

0.77 

South 

Block 

6360 

6364 

6368 

4.52 

4.34 

5.59 

0.84 

0.82 

1.10 

0.92 

0.91 

1.14 

Using  the  RMS  differences  of  Table  7  in  equation  (2-107)  we  obtain 
r0  =  6360.248  km  for  the  NB  and  r0  =  6362.312  km  for  _the  SB  solution. 
The  results  of  the  NB  and  the  SB  solution  with  r0  =  r0  are  shown  in 
Table  8. 
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Table  8.  RMS  Differences  Between  Predicted  and  Control  Values  with 
the  Asymmetric  Kernel  Approach,  Both  Ag  and  ((,tj) 
Observed  and  With  the  Optimal  Radius  of  the  Geosphere, 
by  0S5  Block.  Bjerhammar  Method. 


BLOCK 

RMS  DIFFERENCES  1 

SUB-BLOCK 

1TTCT75H 

mil 

nm 

3.58 

0.77 

nrei 

■db 

4.04 

0.84 

0.58 

2.54 

0.79 

H 

3.47 

0.91 

0.57 

H 

4.32 

0.77 

EEJ 

0.93 

5 

1.26 

1.27 

6 

0.92 

1.04 

7 

4.01 

0.46 

0.34 

8 

2.94 

HE) 

0.73 

From  Table  8  one  can  see  that  using  both  Ag  and  observations 

one  can  predict  Ag  on  the  average  to  about  4  mgals  and  (  and  to 
about  018.  These  results  can  vary  from  2.54  mgals  (sub-block  #2)  to 
7.35  mgals  (sub-block  #6)  for  Ag,  0146  (sub-block  #7)  to  1'126  (Bub-block 
#5)  for  (  and  0134  (sub-block  §7)  to  1*127  (sub-block  #5)  for  tj. 

5.2.2.3  Prediction  Using  Only  Vertical  Deflection  Data 

In  the  case  where  only  vertical  deflections  are  observed  the  least 
squares  solution  applies  with  as  many  degrees  of  freedom  as  observed 
deflection  pairs.  This  type  of  solution  is  given  by  (2-46)  through 
(2-48).  The  elements  of  the  design  matrix  G  are  given  by  (2-32). 
Results  with  three  radii  for  both  the  NB  and  the  SB  solutions  are  given 
in  Table  9. 


Table  9.  RMS  Differences  Between  Predicted  and  Control  Values  with 
the  Asymmetric  Kernel  Approach  and  Only  (( ,ij )  Observed. 
Bjerhammar  Method. 


SOLUTION 

BTC31 

li  icraJI 

ram 

mil 

North 

Block 

6350 

6360 

6365 

1.16 

1.47 

1.84 

ESS 

South 

Block 

6340 

6350 

6360 

58.57 

9.46 

9.19 

1.18 

1.06 

1.64 

1.14 

0.99 

1.19 
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The  RMSs  differences  of  Table  9  in  (2-107)  yield  r0  =  6354.221  km  for 
the  NB  and  r«  =  6350.352  km  for  the  SB  solution.  Using  these  r0  values 
one  gets  the  results  in  Table  10. 


Table  10.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Asymmetric  Kernel  Approach,  Only  ((,17) 
Observed  and  With  the  Optimal  Radius  of  the  Geosphere. 
Bjerhammar  Method. 


BLOCK 

TERENCES  I 

SUB-BLOCK 

Ag(agal) 

tn 

nC) 

North 

13.19 

1.25 

1.06 

1 

11.98 

1.06 

0.91 

2 

8.15 

0.76 

0.74 

3 

16.00 

1.00 

0.80 

4 

12.85 

1.79 

1.50 

South 

9.19 

1.08 

0.99 

5 

10.62 

1.58 

1.25 

6 

7.34 

1.11 

0.89 

7 

10.20 

1.00 

1.35 

8 

7.25 

0.75 

0.85 

From  Table  10  one  can  see  poor  A g  predictions.  On  the  other  hand 
(  was  predicted  to  about  1"17  and  17  to  about  1"03  on  the  average,  with 
variations  from  0175  to  t!79  for  (  and  0174  to  1*150  for  17. 

5.2.3  The  Symmetric  Kernel  Approach 

In  this  series  of  tests  with  the  SK  approach  the  Dirac  Impulses  will 
be  located  at  the  nadir  points  of  the  observations.  Also,  the  optimal 
geosphere  radius  will  be  computed  via  (5-4)  in  order  to  use  it  for  the 
final  solution  as  described  in  Section  5.2.2. 

5.2.3.1  Prediction  Using  Only  Gravity  Anomaly  Data 

In  this  case  the  exact  solution  as  described  by  (2-44)  and  (2-45) 
applies.  The  elements  of  the  design  matrix  G  are  given  by  (2-41).  The 
results  of  both  the  NB  and  the  SB  solutions  with  three  different  radii 
are  shown  in  Table  11. 

Using  the  RMS  differences  of  Table  11  in  equation  (2-107)  one 
obtains  r0  -  6366.339  km  for  the  NB  and  r0  =  6365.267  km  for  the  SB 
solution.  The  results  of  the  NB  and  the  SB  solutions  with  r0  -  r0  are 
shown  in  Table  12. 
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Table  11.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Symmetric  Kernel  Approach  and  Only  Ag 
Observed.  Bjerhammar  Method. 


SOLUTION 

r0(km) 

Ag(mgal) 

tn 

in 

North 

Block 

6360 

6365 

6370 

5.31 

3.62 

3.59 

2.03 

0.95 

1.11 

1.42 

0.96 

1.61 

South 

Block 

6360 

6365 

6370 

4.00 

3.82 

4.78 

1.64 

0.95 

0.99 

1767 

1.13 

1.58 

Table  12.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Symmetric  Kernel  Approach,  Only  Ag  Observed 
*  and  the  Optimal  Radius  of  the  Geosphere  r0.  Bjerhammar 

Method. 


BLOCK 

RMS  DIFFERENCES  1 

SUB-BLOCK 

Ag(agal) 

*»(") 

North 

3.29 

0.92 

0.97 

1 

4.47 

0.88 

0.82 

2 

2.56 

0.66 

1.00 

3 

3.03 

1.32 

0.75 

4 

3.33 

0.80 

1.15 

South 

3.81 

0.94 

1.13 

5 

3.57 

1.36 

1.44 

6 

5.97 

0.96 

1.23 

7 

3.30 

0.58 

0.79 

8 

2.44 

0.78 

0.92 

In  this  case  one  can  see  in  Table  12  that  the  good  sub-block  118  in 
the  sense  of  Table  3  gave  the  best  Ag  predictions  whereas  the  best  ( 
and  i)  predictions  were  performed  at  sub-block  117.  Overall,  with  only 
Ag  observations  the  SK  approach  on  the  average  predicted  Ag  to  about 
3.6  mgals,  (  to  about  019  and  tj  about  l'!05. 

S.2.3.2  Prediction  Using  Both  Gravity  and  Vertical  Deflection  Data 

In  the  event  that  both  Ag  and  ((,*))  are  observed  the  least-squares 
solution  given  by  (2-46)  to  (2-48)  applies,  and  one  has  as  many  degrees 
of  freedom  as  observed  deflection  pairs.  The  elements  of  the  design 
matrix  G  are  given  by  (2-41)  for  Ag  observations  and  by  (2-42)  for  ((,rj) 
observations.  Table  13  shows  the  results  for  this  case  with  three 
different  radii. 
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Table  13.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Symmetric  Kernel  Approach  and  Both  Ag  and 
Observed.  Bjerhammar  Method. 


SOLUTION 

ro(km) 

f(") 

v(") 

North 

Block 

6363 

6365 

6367 

3729  s 

3.38 

3.67 

0.86 

0.81 

0.72 

0.72 

0.72 

0.70 

South 

Block 

6365 

6367 

6369 

4.31 

4.55 

4.20 

0.88 

0.84 

0.82 

0.89 

0.93 

0.91 

The  ^application  of  (2-107)  and  (5-4)  with  the  results  of  Table  13 
yielded  r0  ~  6362.867  km  for  the  NB  and  r0  =  6368.049  km  for  the  SB 
solution.  The  results  of  the  NB  and  SB  solutions  with  r0  =  r0  are 
shown  in  Table  14. 


Table  14.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Symmetric  Kernel,  Both  Ag  and_  Observed 

and  the  Optimal  Radius  of  the  Geosphere  r0.  Bjerhammar 
Method. 


BLOCK 

RMS  DIF! 

TsRENCES  I 

SUB-BLOCK 

Ag(mgal) 

(H 

vD 

North 

3.26 

0.87 

0.72 

1 

3.90 

0.67 

0.54 

2 

2.49 

0.77 

0.71 

3 

3.25 

1.15 

0.83 

4 

3.45 

0.85 

0.75 

South 

4.33 

0.82 

0.91 

5 

3.58 

1.23 

1.23 

6 

7.01 

0.90 

1.01 

7 

3.83 

0.48 

0.37 

8 

2.85 

0.52 

0.72 

From  Table  14  one  observes  Ag  to  be  predicted  to  about  3.5  mgals,  ( 
to  about  0185  and  ij  to  about  018  on  the  average. 

5.2.3.3  Prediction  Using  Only  Vertical  Deflection  Data 

When  only  vertical  deflections  are  observed  the  least-squares 
solution  applies  with  as  many  degrees  of  freedom  as  observed 
pairs.  This  type  of  solution  is  given  by  (2-46)  to  (2-48)  and  the 
elements  of  the  design  matrix  G  are  given  by  (2-42).  Both  the  NB  and 
the  SB  solutions  with  three  different  radii  are  shown  in  Table  15. 
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Table  IS.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Symmetric  Kernel  Approach  and  Only  ((>17) 
Observed.  Bjerhammar  Method. 


SOLUTION 

rQ(ka) 

«") 

vC) 

North 

Block 

6355 

6360 

6365 

50.88 

22.98 

10.19 

1.04 

1.12 

1.36 

1.03 

1.00 

1.18 

South 

Block 

"3533" 

6360 

6365 

“503 

13.05 

8.38 

of 

0.97 

1.48 

03T 

0.99 

1.12 

Using  the  RMS  differences  of  Table  15  in  (2-107)  yields  r0  = 
6359.982  km  for  the  NB  and  r*>  -  6361.017  km  for  the  SB  solution.-  The 
NB  and  SB  solutions  with  these  optimal  radii  are  shown  in  Table  16. 


Table  16.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Symmetric  Kernel,  Only  ((ty)  Observed  and  r0  - 
r0.  Bjerhammar  Method. 


BLOCK 

RMS  DIFFERENCES  I 

SUB-BLOCK 

Ag(agal) 

tn 

vn 

North 

23.05 

1.12 

1.00 

1 

20.82 

0.97 

1.08 

2 

17.16 

0.79 

0.74 

3 

27.34 

0.99 

0.77 

4 

21.65 

1.51 

1.25 

South 

9.84 

1.05 

0.99 

5 

11.35 

1.53 

1.22 

6 

8.25 

1.08 

0.  ifO 

7 

10.69 

1.01 

1.33 

8 

7.93 

0.73 

0.86 

From  Table  16  one  can  see  that  Ag  were  poorly  predicted  with  only 
observations.  Also,  the  best  £  predictions  were  performed  at 
sub-block  118  whereas  sub-block  #2  yielded  the  best  1 j  predictions. 

5.2.4  Comments  on  the  Results  of  the  Asymmetric  and  Symmetric  Kernel 
Approaches 

Comparison  of  Tables  6  and  8  shows  that  in  the  AK  approach  when 
observations  are  introduced  vertical  deflection  predictions  are 
improved  by  about  013,  whereas  the  gravity  anomaly  predictions  were 
downgraded  by  about  0.6  mgals.  Furthermore,  inspection  of  Tables  12 
and  14  yields  very  similar  comparison  for  the  SK  approach  for  the 
introduction  of  ((,»})  observations. 


Comparison  of  Tables  8  and  10  for  the  AK  and  14  and  16  for  the  SK 
approach  demonstrates  deterioration  of  both  the  Ag  and  the  ((,17) 
predictions  when  no  Ag  observations  are  used. 

Comparison  of  Tables  6  to  12,  8  to  14  and  10  to  16  shows  that  the 
results  of  the  predictions  with  the  AK  and  the  SK  approaches  are 
practically  identical.  Actually!  the  only  difference  in  the  two  approaches 
is  the  radius  that  yielded  the  optimal  results.  Table  17  shows  the 
optimal  radii  for  the  two  approaches. 


Table  17.  Optimal  Radii  in  km  for  the  AK  and  the  SK  Approaches  in 
Both  the  NB  and  the  SB  Solutions  with  Different 
Observation  Types.  Bjerhammar  Method. 


Type  of 
Observations 

Optimal  rad] 
Asysmetn 

.us  r0  (km) 
lc  Kernel 

Optimal  radius  r0  (km) 
Symmetric  Kernel 

North  Block 

South  Block 

North  Block 

South  Block 

Ag 

Ag  and  (f.tj) 

_ LLaL _ 

■ 

6361.562 

6362.312 

6350.352 

6366.339 

6362.867 

6359.982 

6365.267 

6368.049 

6361.017 

Prom  Table  17  one  can  see  that  larger  radii  yield  the  optimal  results 
in  the  SK  than  in  the  AK  approach. 


5.2.5  Dirac  Impulses  on  a  Grid 

Up  to  this  point  the  Dirac  Impulses  were  located  at  the  nadir  points 
of  the  observations.  However,  one  potential  location  for  the  Impulses  is 
on  a  grid  at  the  surface  of  the  geoBphere.  In  this  case  it  holds  that  r  t 
in  (2-38)  is  equal  to  r0  and  therefore  t  is  the  same  for  both  the  AK  and 
the  SK  approach.  Consequently  their  respective  formulae  namely  (2-28) 
and  (2-41)  for  Ag  and  (2-32)  and  (2-42)  for  (f ,ij)  become  identical.  This 
scheme  of  computing  Ag*  on  a  grid  was  tested  for  four  different  grid 
sizes  for  both  the  NB  and  the  SB.  The  grids  were  selected  with  two 
considerations  in  mind.  The  first  one  was  to  have  integer  minutes  in 
the  mesh  size.  The  second  one  was  to  have  less  number  of  grid 
vertices  (unknowns)  than  observations  so  that  the  least  squares  solution 
as  given  by  (2-46)  to  (2-48)  be  applicable.  Also,  it  should  be  kept  in 
mind  that  very  coarse  grids  are  not  desirable  since  information  that 
exists  on  the  data  cannot  be  transferred  to  Ag*  values  and  the  resulting 
predictions  become  inaccurate.  The  selected  grids  are  shown  in  Table 
18. 
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Table  18.  Details  of  the  Four  Grids  Used  at  the  White  Sands  TeBt 
Area. 


GRID 

GRID  CELL  SIZE 

«  OF  VERTICES 
IN  4  DIRECTION 

#  OF  VERTICES 
IN  A  DIRECTION 

TOTAL  # 

OF  VERTICES 

1 

6'x4' 

15 

22 

2 

6'x6' 

15 

15 

3 

7'x7' 

13 

13 

4 

12'xl2' 

8 

8 

64 _ I 

The  grids  of  Table  18  were  used  in  two  cases.  One  with  only  Ag 
observations  and  one  with  both  Ag  and  observations.  In  the  case 

of  only  (f  ,tj)  observations  two  grids  were  used.  The  criteria  were  the 
same  as  the  ones  in  the  selection  of  the  grids  of  Table  18.  The  first 
grid  had  a  grid  cell  size  of  7'xl2'  with  13  and  8  vertices  in  the  latitude 
and  longitude  directions  respectively  and  a  total  of  104  vertices.  The 
second  grid  was  identical  to  grid  *4  of  Table  18.  The  computational 

scheme  will  be  to  use  three  arbitrary  values  for  r0  and  record  the 

resulting  RMS  differences  of  control  minus  predicted  values  for  Ag,  { 

and  if.  Subsequently,  these  values  will  be  used  in  conjunction  with 

(2-107)  and  (5-2)  to  yield  the  optimal  geosphere  radius.  The  optimal 
radius  of  the  geosphere  is  used  in  the  final  solution. 

S.2.S.1  Prediction  Using  Only  Gravity  Anomaly  Data 

In  this  case,  the  elements  of  the  design  matrix  G  are  given  by 
either  (2-28)  or  (2-41).  The  results  of  both  the  NB  and  the  SB  solutions 
are  shown  in  Table  19. 


t  Table  19.  RMS  Differences  Between  Predicted  and  Control  Values 

with  the  Four  Grids  and  Only  Ag  Observed.  Bjerhammar 
Method. 


Optimal 
radius (ka) 

NORTH  BL0C1 

S 

m 

immcSl 

S0UT1 

I  BLOCK  i 

OTCTBl 

nsm 

Rim 

rmm 

rmw 

6347.269 

8.86 

6.07 

10.69 

mm 

6353.852 

— 

051 

1.201 

3.21 

1.55 

1.69 

m 

6351.250 

4.39 

0.87 

3.41 

1.24 

U 

6351.250 

4.82 

1.08 

1  6330.576 

4.85 

K BuJ 

1.01 

mm 

6347.500 

6.49 

EEa 

1.49 

From  Table  19  one  can  see  that  gravity  anomaly  predictions  are  best 
performed  with  grid  12.  Furthermore,  the  best  ((,»|)  predictions  were 
performed  with  grid  #4  for  the  NB  and  with  grid  13  for  the  SB.  Most 
importantly,  with  gravity  data  alone,  the  downward  continuation  on  a 
grid  can  yield  similar  results  to  the  ones  obtained  with  the  downward 
continuation  to  the  nadir  points  of  the  observations  (compare  with 
results  of  Tables  6  and  12). 


i 
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S.2.5.2  Prediction  Using  Both  Gravity  and  Vertical  Deflection  Data 

In  this  case  the  elements  of  the  design  matrix  G  are  given  by  either 
(2-28)  and  (2-32)  or  (2-41)  and  (2-42).  The  results  of  both  the  NB  and 
the  SB  solutions  with  the  four  grids  of  Table  18  are  shown  in  Table  20. 


Table  20.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Four  Grids  and  Both  Ag  and  (f  ,tj)  Observed. 
Bjerhammar  Method. 


Optimal 
radius (km) 

NORTH  BLOC] 

I 

soun 

9  BLOCK  1 

mrrm 

R0I 

now 

ETW73 U 

(101 

CT01 

6356.389 

KID 

0.75 

mm 

6350.143 

3.57 

BED 

0.73 

6349.166 

wlr  la 

0.68 

0.61 

E9 

6341.548 

3.94 

0.63 

0.86 

6346.167 

0.65 

mm 

6345.833 

4.69 

0.62 

0.92 

6326.624 

0.72 

n 

6325.686 

6.28 

1.29 

From  Table  20  one  observes  that  grid  #2  yielded  the  best  Ag 
predictions  in  the  NB  whereas  grid  il  was  the  favorite  for  the  SB.  As 
far  as  the  meridional  deflection  predictions  are  concerned  grid  #3  gave 
the  best  results.  However,  in  terms  of  tj  predictions  grid  tt2  performed 
best  in  the  NB  whereas  grid  #1  performed  best  in  the  SB.  Comparing 
Tables  19  and  20  one  sees  that  the  introduction  of  vertical  deflection 
observations  resulted  in  improved  predictions  in  all  cases.  Finally, 
comparison  of  Table  20  to  Tables  8  and  14  reveals  slightly  better  results 
from  the  downward  continuation  onto  a  grid.  However,  the  downward 
continuation  onto  a  grid  has  the  drawback  of  having  to  try  different 
mesh-sizes  in  order  to  get  the  best  predictions,  which  is  impossible  in 
the  absence  of  control  data. 

5.2.5.3  Prediction  Using  Only  Vertical  Deflection  Data 

For  this  application  the  elements  of  the  design  matrix  G  are  given 
by  either  (2-32)  or  (2-42).  The  results  of  the  NB  and  the  SB  solutions 
are  shown  in  Table  21. 


Table  21.  RMS  Differences  Between  Predicted  and  Control  Values 
with  Two  GridB  and  Only  (f,»j)  Observed.  Bjerhammar 
Method. 


iESSmSI 

NORTH  BLOCK 

GRID 

SOUTH  BLOCK  1 

ESi?ig51 

BTOSJ1 

H0i 

R01 

# 

EEJSI 

Ag(mgal) 

Li") 

*|C) 

6330.232 

5424.63 

123.57 

41.82 

1 

6347.498 

10708. 15 

293.05 

66.84 

6347.420 

314.77 

4.98 

3.66 

2 

6348.443 

332.27 

7.21 

4.33 
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From  Table  21  one  observes  very  poor  predictions  for  both  grids 
with  only  (f,ij)  observations. 


5.2.6  The  Best  Ag  and  (f.w)  Predictions 

The  best  gravity  anomaly  predictions  were  obtained  with  the 
Asymmetric  Kernel  (Table  6)  approach  and  with  only  gravity  anomalies  as 
observations.  The  Dirac  Impulses  were  located  on  the  geosphere  at  the 
nadir  points  of  the  observations. 

Inspecting  Table  6  in  light  of  the  representation  factors  of  Table  3 
we  see  that  even  though  the  SB  solution  was  expected  to  be  better  than 
the  NB  one  this  was  not  the  case.  Actually  they  turned  out  about  the 
same.  Furthermore,  from  Table  3  one  sees  that  sub-block  #6  should 
yield  very  good  predictions,  which  was  not  the  case  as  Table  6  shows. 
Also,  even  though  sub-block  #2  is  characterized  poor  in  Table  3,  it 
yielded  good  predictions.  These  somewhat  conflicting  results  force  one 
to  look  at  the  individual  results  at  each  station. 

Figures  26,  27,  28  and  29  show  the  differences  control  minus 
predicted  value  for  each  gravity  control  station  at  the  four  CT.5xO‘.5 
sub-blocks  of  the  NB  and  the  SB  solution  respectively  in  the 
background  of  the  gravity  data.  These  differences  are  from  the  AK 
approach  with  only  gravity  observations  and  the  Dirac  Impulses  located 
at  the  nadir  points  of  the  observations.  In  these  Figures  the  gravity 
control  stations  are  indicated  by  x  and  the  convention  is  that  a  bar 
above  the  x  indicates  a  positive  difference  whereas  a  bar  below  the  x 
indicates  a  negative  difference. 
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Figure  28.  Comparison  by  Station  of  Control  and  Bjerhammar-Predicted 
Gravity  Anomalies  at  the  South  Block  of  the  Teat  Area. 
Sub-blocks  #5  and  6. 
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Figure  29.  Comparison  by  Station  of  Control  and  Bjerhammar-Predicted 
Gravity  Anomalies  at  the  South  Block  of  the  Test  Area. 
Sub-blocks  17  and  8. 
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Figures  26,  27,  28  and  29  indicate  that  the  majority  of  the  large 
differences  occur  at  areas  with  insufficient  data  coverage.  However, 
there  are  exceptions.  For  example,  in  block  #2,  station  117778  has  a 
difference  of  -5.85  mgals  whereas  for  station  #7777  this  difference  is 
-0.08  mgals  and  these  stations  are  only  about  700  m  apart.  Similarly,  in 
block  #3,  station  #6717  which  is  between  stations  6716  and  6718  and  only 
a  few  km  away  from  either  one  has  a  differences  of  5.35  mgals  whereas 
the  other  two  stations  have  differences  of  0.22  mgals  and  -0.32  mgals 
respectively.  Also,  Figures  26,  27,  28  and  29  demonstrate  that  the 
method  can  operate  fairly  well  in  clusters  provided  sufficient  data 
coverage  is  present. 

As  far  as  the  vertical  deflection  predictions  are  concerned,  the  best 
results  were  obtained  when  both  gravity  and  deflection  observations 
were  included.  The  Dirac  Impulses  were  located  on  grid  #3  on  the 
geosphere  (Table  20).  The  results  of  this  solution,  by  (0*.5x0’.5) 
sub-block  are  given  in  Table  22. 

Table  22.  RMS  Differences  Between  Predicted  with  Ag  and  Uji)) 
Observed  and  the  Optimal  Radius  of  the  Geosphere  r0. 
Dirac  Impulses  on  Grid  #3.  Bjerhammar  Method. 


BLOCK 

RMS  DIFFERENCES 

SOB-BLOCK 

ETT035I1 

ran 

mn 

3.56 

0.65 

0.65 

aBH 

4.07 

0.66 

0.54 

3.16 

0.63 

K9 

3.03 

0.83 

0.44 

Wi 

4.34 

0.49 

f»Hr£f 

4.69 

0.62 

0.92 

■m 

0.57 

0.89 

fsa 

4.06 

0.70 

0.80 

nmm 

6.25 

0.59 

1.03 

MBjHk 

2.90 

0.51 

1.04 

The  results  of  Table  22  in  light  of  Table  3  are  conflicting  in  this 
case  also.  For  example  the  "Medium”  sub-block  #4  yielded  the  best 
meridional  deflection  predictions  and  the  prime  vertical  deflection 
predictions  of  the  "Good"  sub-block  #8  were  the  worst  »j  predictions. 
However,  the  "Good"  sub-block  #8  yielded  the  best  gravity  anomaly 
predictions  and  the  second  best  meridional  deflection  predictions. 
Figures  30,  31,  32  and  33  show  the  differences  control  minus  predicted 
value  for  each  vertical  deflection  control  station  at  the  four  O'. 5x0*. 5 
aub-blocks  of  the  NB  and  the  SB  solution  in  the  background  of  the 
gravity  observations.  These  differences  are  from  the  solution  where 
both  Ag  and  were  observed  and  the  Dirac  Impulses  were  located  on 

grid  #3  at  the  surface  of  the  geosphere.  In  these  Figures  the  vertical 
deflection  control  stations  are  indicated  by  x  and  the  convention  is  that 
a  bar  above  the  x  indicates  a  positive  difference  in  (,  a  bar  to  the 
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Figure  31.  Comparison  by  Station  of  Control  and  Bjerhammar-Predicted 
Deflections  of  the  Vertical  of  the  North  Block  of  the  Test 
Area.  Sub-blocks  83  and  4. 
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Figure  33.  Comparison  by  Station  of  Control  and  Bjerhammar-Predicted 
Deflections  of  the  Vertical  at  the  South  Block  of  the  Test 


Area.  Sub-blocks  #7  and  8. 
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right  indicates  a  positive  tj  difference  whereas  a  bar  below  x  indicates 
a  negative  difference  in  (  and  a  bar  to  the  left  indicates  a  negative 
difference  in  tj. 

Figures  30  through  33  illustrate  that  the  majority  of  the  large 
discrepancies  between  control  and  predicted  values  occur  at  areas  poor 
in  data  coverage  and  rich  in  terrain  variations.  For  example  the  hilly 
to  flat  sub-block  #8  with  good  Ag  and  (f,7?)  observation  coverage 
yielded  the  best  Ag  and  second  best  (  predictions.  Moreover,  in  the 
"poor"  sub-block  U2,  station  358  had  (  predicted  within  011  and  -tj 
within  012  due  to  the  presence  of  gravity  observation  stations  7800, 
7801  and  7802. 


5.2.7  Errors  of  Predictions 


In  the  solution  with  optimal  radius  of  the  geosphere  for  each 
variation  of  the  predictor  the  standard  deviations  of  the  predicted 
values  were  computed  according  to  (2-57).  A  close  inspection  of  these 
standard  deviations  indicates  that  they  cannot  be  considered  a  safe 
indicator  of  the  quality  of  the  predictions.  This  is  to  say  that  many 
poorly  predicted  quantities  are  associated  with  small  standard 
deviations  and  many  very  well  predicted  quantities  are  associated  with 
large  standard  deviations. 


5.2.8  Conclusions  from  the  Bjerhammar  Predictor 

At  first  the  optimal  radius  r0  of  the  geosphere  could  not  be 
computed  from  the  data.  The  results  of  two  methods  to  perform  this 
computations  were  discouraging.  Therefore  the  sa  method  was  used  at 
which  the  RMS  discrepancies  (control  minus  predicted)  in  Ag,  (  and  rj 
were  considered  a  second  order  polynomial  in  r0. 

It  gravity  predictions  are  required,  then  use  only  gravity  data  and 
position  the  Dirac  Impulses  at  the  nadir  points  of  the  observations.  It 
is  not  significant  whether  the  Asymmetric  or  the  Symmetric  Kernel  is 
used  in  terms  of  the  quality  of  the  predictions.  The  only  requirement 
to  get  the  same  prediction  quality  from  both  the  AK  and  the  SK  is  to 
associate  the  AK  with  radii  about  6362  km  and  the  SK  with  radii  about 
6366  km.  The  exact  value  for  the  optimal  r0  should  be  dictated  by  the 
specific  data  set  with  the  s*  method.  In  the  event  that  no  control 
data  exist  in  an  area,  the  observations  can  be  separated  in  two  groups, 
one  of  which  will  play  the  role  of  observations  and  the  other  one  the 
role  of  control  data  so  that  an  optimal  r„  can  be  computed  by  the 
8a-method. 

If  vertical  deflection  predictions  are  sought,  then  include  ((,tj) 
observations  with  the  gravity  data  and  place  the  Dirac  Impulses  on  a 
grid.  When  performing  computations  on  scalar  computers  select  the 
grid  cell  size  keeping  in  mind  that  finer  grids  are  more  CPU  time 
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consuming  withoug  being  necessarily  more  accurate.  In  the  event  that 
a  super-computer  is  available,  the  grid  cell  size  may  not  introduce  a 
problem  in  terms  of  CPU  time.  This  is  the  case  with  the  Cray  X-MP/24 
which  was  used  for  our  solutions.  As  for  the  grid  size,  the  White 
Sands  Test  Area  seems  to  indicate  that  it  can  be  about  twice  the 
angular  distance  between  gravity  observations.  A  radius  of  about  6346 
km  yielded  optimal  results. 

As  far  as  the  R  factor  (eq.(5-l>)  is  concerned,  one  may  conclude 
that  it  is  of  limited  importance.  For  example,  from  Table  22  one  can 
see  that  the  "Medium"  sub-block  #4  (see  Table  3)  yielded  the  best  f 
predictions,  the  "Good"  sub-block  #8  yielded  the  worst  i)  predictions, 
even  though  the  "Good"  sub-block  #8  yielded  the  best  Ag  predictions 
and  the  "Poor"  sub-block  #7  yielded  very  poor  t>  predictions. 

For  terrain  height  variations  from  the  mean  height  of  70  m  to  220 
m,  gravity  data  density  of  20  to  30  km2  per  station  and  standard 
deviations  of  the  data  in  the  order  of  2  mgals,  the  method  can  predict 
Ag  within  4  mgals  and  within  1".  If  vertical  deflection 

observations  as  dense  as  140  to  240  km*  per  station  and  as  accurate  as 
013  are  added,  then  vertical  deflections  can  be  predicted  to  O'. 7  or 
better. 


5.3  Results  of  the  Hardy  Method 

5.3.1  Tests  of  Optimal  Geosphere  Radius  Computation 

The  radius  of  the  internal  sphere  is  as  important  with  Ilardy’s 
predictor  as  it  is  with  Bjerhammar’s  predictor.  The  coupling  effect  of 
r0  mentioned  in  subsection  5.2.1  is  present  here  also.  Therefore  an 
optimal  value  for  the  radius  was  attempted  to  be  computed  from  the 
data  with  this  predictor  as  well. 

The  same  two  methods  tested  with  the  Bjerhammar  predictor  were 
tested  with  the  Hardy  predictor  also.  The  first  one  was  to  use  384  Ag 
observations  in  the  NB  and  solve  for  384  biharmonic  sources  c  (  plus 
the  radius  r0  of  the  geosphere.  An  approximate  value  of  6350  for  r0 
resulted  after  4  iterations  in  an  adjusted  value  of  6381585.4  *  11.4  m. 
The  residuals  were  in  the  order  of  10“4  mgals;  VTPV  was  10~*  and  the 
standard  deviations  of  the  biharmonic  sources  cf  were  larger  than  the 
c(  values.  An  approximate  value  of  6360  km  for  r0  resulted  after  10 
iterations  in  an  adjusted  value  of  6359324.1  *  1018.1  m.  The  residuals 
were  of  the  order  of  10-3  mgals,  VTPV  was  10-4  and  the  standard 
deviations  of  ct  were  larger  than  the  cf  values.  The  order  of 
magnitude  of  the  residuals  and  of  VTPV  can  be  explained  by  the 
absence  of  redundant  observations.  The  large  standard  deviations  of 
the  C|  values  as  well  as  the  standard  deviation  of  the  adjusted  ru 
simply  stress  that  the  resulted  adjusted  values  have  not  been 
accurately  determined.  From  the  above  results  one  cannot  conclude  in 
favor  of  a  meaningful  r0  computation. 


The  second  attempt  was  the  data  separation  method  (equation 
(2-80)  through  (2-102) ).  The  method  was  tested  in  the  NB  and  the 
results  of  the  solution  (iteration  #0)  as  well  as  the  six  iterations 
required  for  convergence  are  shown  in  Table  23.  In  Table  23,  N  is  the 
normal  matrix  (of  dimensions  (lxl))  and  U  is  the  right  hand  side  vector 
(of  dimensions  (lxl))  of  the  normal  equations. 


Table  23.  Data  Separation  Method  of  Computing  rD.  Hardy  Method. 


:  14^*2 

BrTi 

RIRI 

N 

U 

6r0(km) 

HVBI 

mm 

6360.000 

0.000  006  657 

0.009  826  926 

-1.476 

6358.524 

mm 

6358.524 

0.000  005  277 

-0.011  218  644 

2.126 

6360.650 

Sip 

6360.650 

0.000  004  396 

-0.024  550  281 

5.585 

6366.235 

mm 

6366.235 

0.000  027  438 

-0.148  932  793 

5.428 

6371.663 

SB 

6371.663 

0.085  465  724 

-38.874  943  956 

0.455 

6372.118 

KB 

6372.118 

1.708  862  581 

-160.625  401  236 

0.094 

6372.212 

m 

6372.212 

415.993  395  432 

-106.240  406  225 

0.000 

6372.212 

The  adjustment a yielded  rj  =  6372.212  *  5xl0-5  km.  Also,  it  yielded 
VrPV  =  5xl04  and  <r„  s  11.7.  From  Table  23  one  observes  that,  after 
the  second  iteration,  both  N  and  U  are  increasing  in  absolute  value. 
However,  the  correction  6 r0  tends  to  zero  after  the  third  iteration. 
The  adjusted  value  of  the  internal  sphere  was  greater  than  the  mean 
Earth  radius  of  6371  km.  Using  r0  =  rg  =  6372.212  km  resulted  in  RMS 
discrepancies  control  minus  predicted  of  17.37  mgals  for  Ag,  93(11  for  ( 
and  323181  for  rj.  Conclusively,  neither  one  of  the  two  methods  appears 
to  be  able  to  compute  an  optimal  r0  value.  As  a  result,  the  s2  method 
should  be  used  for  r0  computations  with  the  Hardy  predictor. 

5.3.2  Biharmonic  Sources  at  the  Nadir  Points  of  the  Observations 

In  this  series  of  tests  the  biharmonic  sources  cs  will  be  located  at 
the  nadir  points  of  the  observations.  The  differences  in  the  following 
tables  will  be  control  minus  predicted  and  the  s2  method  will  be  used 
for  optimal  r0  computations.  The  final  solutions  for  each  variation  of 
the  method  will  be  performed  with  the  r0  value  computed  from  (5-4) 
based  on  the  RMS  differences  resulting  from  solutions  with  three 
different  radii. 

5.3.2.1  Prediction  Using  Only  Gravity  Anomaly  Data 

In  the  case  where  only  Ag  are  considered  as  observations  the  exact 
solution  applies  as  given  by  (2-44)  and  (2-45).  The  elements  of  the 
design  matrix  G  are  given  by  (3-43).  The  results  from  both  the  NB 
and  the  SB  solutions  with  three  different  radii  are  shown  in  Table  24. 
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Table  24.  RMS  Differences  Between  Predicted  and  Control  Values 
with  Only  Ag  Observed.  Hardy  Method. 


SOLUTION 

r0(k*> 

fH 

iH 

North 

Block 

6355 

6360 

6365 

4.39 

3.58 

3.03 

8.49 

4.72 

5.41 

6.59 

7.24 

13.37 

South 

Block 

6355 

6360 

6365 

3.90 

3.82 

3.75 

2.57 

2.57 

4.53 

4.43 

7.51 

15.55 

Application  of  (5-4)  with  r$9,  r£,  r$  as  computed  through  (2-107) 
and  the  _RMS  differences  of  Table  24  yielded  r0  =  6363.903  km  for  the 
NB  and  r0  =  6355.948  km  for  the  SB  solution.  The  NB  and  SB  solutions 
with  r0  =  r0  are  shown  in  Table  25. 


Table  25.  RMS  Differences  Between  Predicted  and  Control  Values 
with  Only  Ag  Observed  and  the  Optimal  Radius  of  the 
Geosphere  r0>  Hardy  Method. 


BLOCS 

RMS  DIFFERENCES  ! 

SUB-BLOCK 

Ag(ngal) 

tn 

*in 

North 

3.13 

4.93 

11.36 

1 

4.32 

2.86 

11.42 

2 

2.56 

6.17 

3.95 

3 

2.88 

2.32 

18.44 

4 

2.98 

6.00 

9.52 

South 

3.88 

2.43 

4.77 

5 

3.73 

3.93 

7.24 

6 

5.69 

1.49 

4.97 

7 

3.57 

1.04 

4.18 

8 

2.63 

2.91 

3.12 

Prom  Table  25  we  can  see  that  the  NB  solution  is  slightly  better 
than  the  SB  one  in  terms  of  gravity  predictions.  The  best  Ag 
predictions  were  performed  in  sub-block  #2.  However,  the  SB  solution 
is  better  than  the  NB  one  in  terms  of  vertical  deflection  predictions. 
The  best  f  predictions  were  performed  in  sub-block  *7  and  the  best  ij 
ones  were  performed  in  sub- block  #8. 

5. 3. 2. 2  Prediction  Using  Both  Gravity  and  Vertical  Deflection  Data 

In  the  event  that  both  Ag  and  ({,»»)  are  utilized  as  observed 
quantities  the  least-squares  solution  applies  as  given  by  (2-46) 
through  (2-48)  and  we  have  as  many  degrees  of  freedom  as  observed 
deflection  pairs.  The  elements  of  the  design  matrix  G  are  given  by 
(3-43)  for  Ag  observations  and  by  (3-44)  for  ((,?))  observations.  The 
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solutions  for  both  the  NB  and  the  SB  with  three  different  radii  are 
shown  in  Table  26. 


Table  26.  RMS  Differences  Between  Predicted  and  Control  Values 
with  Both  Ag  and  ($,»j)  Observed.  Hardy  Method. 


SOLUTION 

rn(kn) 

Aft(mgal) 

f(") 

North 

Block 

6360 

6365 

6370 

6.88 

12.73 

20.09 

2.23 

2.43 

3.15 

2.31 

3.09 

3.74 

South 

Block 

6361 

6365 

6368 

7.57 

13.09 

18.09 

2.29 

4.12 

4.83 

2.43 

4.17 

4.92 

Using  the  results  of  Table  26  with  the  s2  method  (equation  (2-107)) 
and  the  RMS  differences  we  get  r0  =  6365.402  km  for  the  NB  and  the 
r0  -  6362.044  km  for  the  SB  solution.  The  r0  value  for  SB  resulted  in 
a  normal  matrix  with  numerically  linearly  dependent  columns  and 
therefore  non-invertible.  Alternatively  the  value  of  6361  km  was 
considered  optimal  for  the  SB.  The  results  from  the  NB  solution  with 
r0  =  rQ  and  the  SB  solution  for  r0  =  6361  km  are  shown  in  Table  27. 


Table  27.  RMS  Differences  Between  Predicted  and  Control  Values 
with  Both  Ag  .and  Observed,  r0  =  6362.186  km  for 

the  NB  and  r0  6361  km  for  the  SB  solution.  Hardy 
Method. 


BLOCK 

RMS  DIF1 

TERENCES  ! 

SUB-BLOCK 

Ag(sgal) 

?n 

7C) 

North 

13.37 

2.45 

3.19 

1 

9.12 

1.65 

1.79 

2 

5.08 

2.32 

3.10 

3 

18.55 

1.84 

2.43 

4 

10.85 

3.29 

4.32 

South 

10.13 

3.70 

3.11 

5 

8.00 

2.76 

1.91 

6 

13.84 

4.85 

4.30 

7 

10.25 

1.04 

0.75 

8 

8.73 

2.69 

1.86 

From  Table  27  we  can  see  that  Ag  can  be  predicted  from  5.08  mgals 
(sub-block  #1)  to  18.55  mgals  (sub-block  #3).  Also,  the  RMS 
•  discrepancies  vary  from  l"04  (sub-block  17)  to  C85  (sub-block  46)  for 
(  and  from  0175  (sub-block  47)  to  4132  (sub-block  44)  for  »j.  Not 
surprisingly,  comparison  of  Tables  25  and  27  reveals  that  the 
introduction  of  vertical  deflection  observations  degrades  the  gravity 
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predictions  whereas  it  improves  vertical  deflection  predictions. 

5.3.2.3  Prediction  Using  Only  Vertical  Deflection  Data 

In  the  case  where  only  vertical  deflections  are  observed  the 
least-squares  solution  applies  with  as  many  degrees  of  freedom  as 
deflection  pairs.  This  type  of  solution  is  given  by  (2-46)  through 
(2-48).  The  elements  of  the  design  matrix  G  are  given  by  (3-44). 
Results  from  both  the  NB  and  the  SB  solutions  with  three  different 
radii  are  shown  in  Table  28. 


Table  28.  RMS  Differences  Between  Predicted  and  Control  Values 
with  Only  (f,7j)  Observed.  Hardy  Method. 


SOLUTION 

To  (km) 

Ag(mgal) 

en 

’?(”) 

North 

Block 

6355 

6360 

6365 

16.33 

15.15 

15.78 

0.97 

1.01 

1.07 

0.99 

1.01 

1.09 

South 

Block 

6350 

6355 

6360 

29.57 

22.90 

22.03 

0.93 

0.91 

0.99 

1.06 

1.01 

1.05 

The  sa  method  of  equation  (2-107)  with  the  RMS  differences  of 
Table  28  yielded  r0  =  6354.698  km  for  the  NB  and  r„  =^6355.676  km  for 
the  SB  solution.  The  NB  and  SB  solutions  with  r„  =  r0  are  shown  in 
Table  29. 


Table  29.  RMS  Differences  Between  Predicted  and  Control  Values 
with  Only  ($,tj)  Observed  and  the  Optimal  Radius  of  the 
Geosphere  r„.  Hardy  Method. 


BLOCK 

RMS  DIFFERENCES 

SUB-BLOCK 

Ag(mgal) 

*in 

North 

16.48 

0.97 

0.98 

1 

11.68 

0.87 

1.16 

2 

5.93 

0.74 

0.77 

3 

23.62 

0.97 

0.88 

4 

10.43 

1.20 

1.09 

South 

12.61 

0.91 

1.01 

5 

6.49 

1.46 

1.23 

6 

20.09 

0.80 

0.93 

7 

10.90 

0.97 

1.37 

8 

12.46 

0.72 

0.86 

From  Table  29  one  can  see  that  the  best  Ag  predictions  are 
performed  at  sub-block  *2.  Sub-blocks  #2,  6  and  8  did  well  for  ( 


predictions  as  did  sub-block  12,  3,  6  and  8  for  i).  Comparison  of  Table 
29  to  Tables  25  and  27  indicates  furthe~  improvement  of  the  ((,tj) 
predictions  and  deterioration  of  the  Ag  p  .ctions  from  the  removal  of 
Ag  observations. 

5.3.3  Biharmonic  Sources  on  a  Grid 


Up  to  this  point  the  biharmonic  sources  were  located  at  the  nadir 
points  of  the  observations  on  the  internal  sphere  (see  Subsection 
5.3.2).  In  the  following  sequence  of  tests  the  biharmonic  sources  will 
be  placed  on  a  grid  at  the  surface  of  the  geosphere.  The  scheme  for 
selecting  the  grids  was  the  same  as  previously  employed  (Subsection 
5.2.4);  the  four  grids  of  Table  18  were  used.  The  computational  scheme 
will  be  to  use  three  arbitrary  values  for  r0  and  record  the  resulting 
RMS  differences  of  control  minus  predicted  values  for  Ag,  (  and  rj. 
Subsequently,  these  values  will  be  used  in  conjunction  with  (2-107) 
and  (5-2)  to  yield  the  optimal  geosphere  radius.  The  optimal  radius  of 
the  geosphere  is  used  in  the  final  solution. 

5.3.3. 1  Prediction  Using  Only  Gravity  Anomaly  Data 

In  this  case  the  least-squares  solution  applies  as  given  by  (2-46) 
to  (2-48).  The  elements  of  the  design  matrix  G  are  given  by  (3-43). 
The  results  of  the  NB  and  the  SB  solutions  are  shown  in  Table  30. 


Table  30.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Pour  Grids  and  Only  Ag  Observed.  Hardy 
Method. 


NORTH  BLOC] 

l 

»ggga| 

SOUTH  BLOCK 

tneraii 

raa 

EH 

mcraM 

warn 

rmm 

6356.965 

23.91 

u 

6356.973 

4.39 

6354.388 

3.21 

4.21 

8.61 

Ea 

6351.499 

4.39 

2.72 

6351.746 

3.37 

7.17 

5.93 

H 

6362.912 

4.77 

8.07 

6342.176 

4.91 

2.88 

n 

6354.502 

6.51 

3.37 

5.15 

From  Table  30  one  sees  that  grids  #2,  3  in  the  NB  and  grids  II,  2, 
3  in  the  SB  yield  satisfactory  Ag  predictions.  However,  the  vertical 
deflections  were  predicted  poorly  from  only  gravity  anomaly  data. 

5.3.3.2  Prediction  Using  Both  Gravity  and  Vertical  Deflection  Data 

Here  the  least-squares  solution  applies  also.  The  elements  of  the 
design  matrix  G  are  given  by  (3-43)  for  gravity  anomaly  and  by  (3-44) 
for  vertical  deflection  observations.  The  results  of  NB  and  SB 
adjustments  are  shown  in  Table  31. 
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Table  31.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Four  Grids  and  Both  Ag  and  (f  ,tj)  Observed. 
Hardy  Method. 


Optimal 
radius (km) 

N0RT1 

I  BLOC] 

l 

B5SHR551 

S0UT1 

1  BLOCK  i 

OTCT5J1 

Rttl 

Clttl 

mrmi 

run 

nei 

6352.264 

4.85 

1.16 

1.39 

wm 

6350.196 

1.19 

6343.865 

4.36 

0.98 

0.86 

6342.498 

4.53 

Wm 

0.98 

6335.444 

4.50 

0.81 

0.77 

El 

6334.255 

5.34 

IxS 

1.08 

6312.865 

5.18 

0.76 

iff 

6326.940 

6.76 

0.93 

1.59 

From  Table  31  one  observes  that  grid  #2  is  the  preferred  choice 
for  Ag  predictions  at  both  the  MB  and  the  SB.  However,  grid  #3 
yielded  best  (£,i))  for  the  NB  and  best  £  for  the  SB.  The  best  ij 
predictions  at  the  SB  were  performed  using  grid  #2. 

Comparison  of  Tables  30  and  31  indicates  that  the  introduction  of 
(£,ij)  observations  deteriorated  the  Ag  predictions  whereas  it  improved 
the  predictions. 

5.3.3.3  Predictions  Using  Only  Vertical  Deflection  Data 


The  two  grids  of  subsection  5.2.4. 3  were  used  in  this  case.  The 
solution  is  of  the  least-squares  type,  and  the  elements  of  the  design 
matrix  G  are  given  by  (3-44).  The  results  for  both  the  NB  and  the  SB 
adjustments  are  shown  in  Table  32. 


Table  32.  RMS  Differences  Between  Predicted  and  Control  Values 
with  the  Two  Grids  and  Only  (£,t>)  Observed.  Hardy 
Method. 


Optiaml 

NORTH  BLOCK 

^  SOUTH  BLOCK  “1 

radius (km) 

mcmnmsm 

rmm 

■SI 

1  :T«  IB 

mm 

EH 

6177.393 

194.41  4.48 

5.39 

1 

208.57 

3.95 

8.25 

6328.837 

777.24  3.28 

3.78 

2 

6329.203 

406.48 

4.66 

2.05 

From  Table  32  one  sees  that  gravity  anomalies  are  predicted 
unacceptably  with  both  grids.  Furthermore,  the  vertical  deflection 
predictions  are  poor.  Also,  from  Table  32  one  sees  that  the  value  of 
the  optimal  radius  of  the  geosphere  for  grid  #1  is  peculiar. 


Comparison  of  Tables  31  and  32  indicates  that  removal  of  Ag 
observations  deteriorated  both  the  Ag  and  the  predictions. 
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5.3.4  The  Best  A*  and  It.n)  Predictions 

As  far  as  the  gravity  anomaly  predictions  are  concerned  the  best 
results  were  obtained  using  only  gravity  observations  and  locating  the 
biharmonic  sources  at  the  nadir  points  of  the  observations.  This 
solution  yielded  RMS  differences  in  the  order  of  4  mgals  for  gravity 
anomalies  and  is  shown  in  Table  25. 

Inspecting  Table  25,  with  the  representation  factor  R  of  Table  3  in 
mind,  we  see  that  block  *8  performed  well  as  -  expected.  Also,  blocks 
42,  3  and  4  performed  well  even  though  they  were  classified  as  not 
very  good.  The  "promising"  sub-block  46  according  to  Table  3  yielded 
the  worst  results.  Figures  34,  35,  36  and  37  show  the  differences 
control  minus  predicted  value  for  each  gravity  control  station  at  eight 
0’.5x0’.5  aub-blocka.  The  convention  for  positive  and  negative  values 
is  the  one  used  in  Figures  26,  27,  28  and  29. 

The  fact  that  the  terrain  type  and  the  data  coverage  influences 
the  predictions  greatly  is  also  illustrated  in  Figures  34  through  37. 
The  problems  of  stations  7778  and  6717  mentioned  at  the  Bjerhammar 
method  exist  with  the  Hardy  predictor.  The  difference  for  station  47777 
is  also  0.10  mgals  whereas  for  7778  it  is  -5.69  mgals  and  the 
discrepancy  for  6717  is  5.27  mgals  whereas  for  6716  it  is  0.14  and  for 
6718  is  -0.26  mgals. 

In  terms  of  vertical  deflection  predictions,  the  best  results  were 
attained  when  both  Ag  and  (f  ,ij)  were  observed  and  the  biharmonic 
sources  were  located  on  grid  43  at  the  surface  of  the  geosphere  (Table 
31).  These  solutions  yielded  good  Ag  predictions  bIbo  and  are  shown  in 
Table  33  by  (O’. 5x0". 5)  sub-block. 
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Figure  37.  Comparison  by  Station  of  Control  and  Hardy-Predicted 
Gravity  Anomalies  at  the  South  Block  of  the  Test  Area. 
Sub-blocks  87  and  8. 


Table  33.  RMS  Differences  Between  Predicted  and^  Control  Values 
with  Ag  and  Observed  and  r0  =  r0.  Biharmonic 

Sources  on  Grid  V3.  Hardy  Method. 


BLOCK 

RMS  DIFFERENCES  ! 

SUB-BLOCK 

«n 

*?(") 

North 

4.50 

0.81 

0.77 

1 

4.56 

0.69 

0.54 

2 

3.68 

0.62 

0.77 

3 

3.85 

1.00 

0.93 

4 

6.14 

0.90 

0.79 

South 

5.34 

0.69 

1.08 

5 

4.77 

0.76 

1.18 

6 

5.03 

0.75 

0.91 

7 

7.00 

0.69 

1.12 

8 

3.37 

0.55 

1.23 

Prom  Table  33  one  can  see  that  the  "Good"  sub-block  #8  yielded 
the  best  f  predictions  and  the  "Medium"  sub-block  #1  yielded  the  best 
7i  predictions.  Most  importantly  from  Table  33  one  sees  that 
predictions  below  the  1"  mark  can  be  performed  with  the  Hardy 
Method.  Figures  38,  39,  40  and  41  show  the  differences  control  minus 
predicted  value  for  the  eight  0:5x015  sub-blocks.  The  convention  for 
positive  and  negative  differences  is  the  one  used  in  Figures  30,  31,  32 
and  33. 

Figures  38  through  41  demonstrate  that  large  differences  mostly 
occur  at  areas  rich  in  terrain  variations  and  poor  in  data  coverage. 
With  the  Hardy  predictor  one  observes  that  sub-blocks  #8  and  7  yield 
the  best  and  worst  (  predictions  respectively  and  that  sub-blocks  ttl 
and  8  yield  the  best  and  worst  ri  predictions  respectively. 
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Figure  39.  Comparison  by  Station  of  control  and  Hardy-Predicted 
Vertical  Deflections  at  the  North  Block  of  the  Test  Area. 
Sub-blocks  *3  and  4. 
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Figure  41.  Comparison  by  Station  of  Control  and  Hardy-Predicted 
Vertical  Deflections  at  the  South  Block  of  the  Test  Area. 
Sub-blocks  #7  and  8. 
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5.3.5  Errors  of  Predictions 


For  every  variation  of  the  predictor  the  standard  deviations  a  of 
the  predictions  were  computed  at  the  solutions  with  the  optimal  radii 
r0  according  to  (2-57).  These  a  values  are  rather  smooth  and  cannot 
indicate  locations  at  which  predictions  are  good  or  not. 

5.3.6  Conclusions  from  the  Hardy  Predictor 

At  first,  the  results  of  the  two  methods  to  compute  r0  from  the 
data  were  not  promising  with^  the  Hardy  Predictor.  As  a  result  the  s2 
method  was  used  to  compute  r0. 

If  gravity  predictions  are  required,  use  only  gravity  data  and 
place  the  biharmonic  sources  at  the  nadir  points  of  the  observations  on 
the  geosphere.  For  White  Sands,  a  radius  of  about  6356  km  to  6363  km 
appears  to  be  optimal. 

If  vertical  deflection  predictions  are  required  use  both  Ag  and 
((,i;)  observations  and  locate  the  biharmonic  sources  on  a  grid  at  the 
surface  of  the  geosphere.  A  radius  of  about  6335  km  yielded  optimal 
results  for  New  Mexico. 

The  overall  result  of  the  tests  of  the  Hardy  Predictor  seems  to  be 
that  Ag  can  be  predicted  to  about  3  to  4  mgals,  (  and  to  about  O'. 8 
0!9. 

5.3.7  Comparison  of  Bjerhammar  and  Hardy  Predictors 

Theoretically  the  predictors  are  very  different.  They  even  assume 
different  behavior  of  the  disturbing  potential  T.  However,  in  practice 
they  yielded  very  similar  results.  The  best  Ag  and  predictions 

were  performed  with  identical  data  requirements  and  downward 
continuation  scheme  and  yielded  similar  results.  A  minor  difference  is 
the  value  of  the  optimal  geosphere  radius.  The  aforementioned  results 
seem  to  stress  that  an  improvement  in  the  predictions  of  both  Ag  and 
((iv)  even  in  mountainous  areas  may  not  result  from  a  theoretical 
breakthrough  but  from  improved  data  coverage. 


5.4  Prediction  Using  Least-Squares  Collocation 

For  comparison  purposes  a  Least-Squares  Collocation  solution  was 
tested  at  the  New  Mexico  Area.  The  model  used  for  the  disturbing 
potential  covariance  function  was  of  the  form  [Kearsley  et  al.,  1985;  p. 
50] 


K(P,Q) 


m 

i 

n=  i  a  i 


ARc2 _ 

(n-1) (n-2) (n+24) 


f  rb21 n+1 

P"(cOSw™> 


(5-5) 
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where  <yPQ  is  the  spherical  distance  between  points  P  and  Q,  rP  and  r0 
are  geocentric  radial  distances  to  P  and  Q  respectively,  RB  is  the 
radius  to  the  Bjerhammar  sphere  and  RE  is  the  mean  Earth  radius. 
The  following  values  were  used  [Kearsley  et  al.,  1985,  p.  50]. 


Re  =  6371  ka,  RB  =  6369.75  ka 


(5-6) 


variance  CqA9  of  the  residual  gravity  anomalies  as  computed  fr 


1137  point  values  was  [Heiskanen  and  Moritz,  1967;  p.  253] 


CoA9  =  Var(VAa)  =  323.82  rngal*  (5-7) 

where  VAg  is  given  by  (4-10). 

For  each  test  two  solutions  were  performed.  One  for  the  NB  and 
one  for  the  SB.  The  first  test  was  to  predict  Ag  and  ($,ij)  from 
gravity  data  alone.  The  results  of  this  test  for  both  the  NB  and  the 
SB  are  shown  in  Table  34. 


Table  34.  RMS  Differences  Between  Predicted  and  Control  Values 
with  Only  Ag  Observed.  Collocation  Solution. 


BLOCK 

i’ERENCES  I 

SUB-BLOCK 

Ag(agal) 

vD 

North 

2.84 

0.76 

0.92 

1 

4.23 

0.83 

0.79 

2 

2.59 

0.63 

0.91 

3 

2.28 

1.02 

0.75 

4 

2.65 

0.61 

1.10 

South 

3.68 

0.80 

1.03 

5 

3.01 

1.00 

1.25 

6 

6.25 

0.86 

1.16 

7 

3.36 

0.57 

0.79 

8 

1.68 

0.67 

0.79 

From  Table  34  one  sees  RMS  differences  from  1.68  to  6.25  mgals  for 
Ag,  0'57  to  1102  for  (  and  01 75  to  1125  for  tj  resulting  from  different 
data  coverage  and  terrain  type  within  the  eight  sub-blocks.  On  the 
average  Ag  was  predicted  to  about  3  mgalB,  (  to  about  018  and  v  to 
about  HO. 

The  second  test  was  to  predict  Ag  and  ((,n)  from  both  Ag  and 
observations.  The  results  of  this  attempt  for  both  the  NB  and  the  SB 
solution  are  shown  in  Table  35. 


-  -v  ; 


i  *  - 
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Table  35.  RMS  Differences  Between  Predicted  and  Control  Values 
with  both  Ag  and  Observed.  Collocation  Solution. 


BLOCK 

RMS  DIFFERENCES  1 

SUB-BLOCK 

in 

*?n 

North 

!  2793 

0.58 

0.56 

1 

4.08 

0.59 

0.48 

2 

2.40 

0.49 

0.60 

3 

2.56 

0.82 

0.31 

4 

3.00 

0.44 

0.68 

South 

3.44 

0.64 

0.66 

5 

3.37 

0.89 

1.06 

6 

4.61 

0.70 

0.59 

7 

3.70 

0.31 

0.37 

8 

1.82 

0.49 

0.59 

From  Table  35  one  sees  RMS  differences  from  1.82  to  4.61  ragals  for 
Ag,  0131  to  0189  for  {  and  0231  to  1106  for  i)  due  to  the  terrain  type 
and  data  coverage  of  the  sub-blocks.  On  the  average  Ag  was 
predicted  to  about  3.3  mgals  and  (  and  rj  to  about  016. 

The  third  test  was  to  predict  Ag  and  (f,v)  from  vertical  deflection 
observations  alone.  The  results  of  this  test  for  both  the  NB  and  the 
SB  solution  are  shown  in  Table  36. 


Table  36.  RMS  Differences  Between  Predicted  and  Control  Values 
with  only  Observed.  Collocation  Solution. 


BLOCK 

RMS  DIF1 

i^RENCES  1 

SUB-BLOCK 

Ag(agal) 

tn 

i?n 

North 

5.41 

0.67 

0.86 

1 

6.02 

0.65 

0.87 

2 

3.75 

0.52 

0.66 

3 

5.91 

0.73 

0.83 

4 

5.46 

0.76 

1.02 

South 

7.36 

0.72 

0.86 

5 

5.92 

1.31 

1.27 

6 

3.49 

0.58 

0.79 

7 

10.70 

0.53 

0.64 

8 

5.47 

0.60 

0.80 

From  Table  36  one  sees  RMS  differences  from  3.49  to  10.70  mgals 
for  Ag,  0152  to  1131  for  (  and  0164  to  1127  for  tj.  The  reason  for  this 
variation  is  the  terrain  type  and  the  data  coverage  in  the  sub-blocks. 
On  the  average  Ag  was  predicted  to  about  6.5  mgals,  (  to  about  017  and 
if  to  about  019. 
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Comparison  of  Tables  34j  35  and  36  shows  that  the  introduction  of 
vertical  deflection  data  slightly  improved  the  predictions  and 

slightly  deteriorated  the  Ag  predictions.  It  is  noteworthy  that  the  RMS 
difference  of  6.25  mgals  at  sub-block  *6  was  improved  to  4.61  mgals  by 
the  introduction  of  (( ,if )  observations.  The  removal  of  gravity  data 
resulted  in  degradation  of  the  Ag  predictions  by  about  3  mgals  and  a 
slight  degradation  of  the  (f  ,tj)  predictions.  In  conclusion,  the  best  Ag 
predictions  were  obtained  from  Ag  data  alone  (Table  34).  This  solution 
is  shown  by  station  in  Figures  42,  43,  44  and  45.  On  the  other  hand, 
the  best  ((,i»)  predictions  are  obtained  from  both  Ag  and  (£,tj)  data 
(Table  35)  and  this  solution  is  shown  in  Figures  46,  47,  48  and  49  by 
station. 


Cm* (LOCK  *9-33 
X*  WAV.  MW.  COMP.  STMT. 
«■  CAAV.  flHOA.  006.  STM. 


Figure  43.  Comparison  by  Station  of  Control  and  Least  Squares 
Collocation  -  Predicted  Gravity  Anomalies  at  the  North  Block 
of  the  Test  Area.  Sub-blocks  #3  and  4. 
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Figure  44.  Comparison  by  Station  of  Control  and  Least  Squares 
Collocation  -  Predicted  Gravity  Anomalies  at  the  South  Block 
of  the  Test  Area.  Sub-blocks  95  and  6. 
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Figure  46.  Comparison  by  Station  of  Control  and  Least  Squares 
Collocation  -  Predicted  Vertical  Deflections  at  the  North 
Block  of  the  Test  Area.  Sub-blocks  #1  and  2. 
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Figure  47.  Comparison  by  Station  of  Control  and  Least  Squares 
Collocation  -  Predicted  Vertical  Deflections  at  the  North 
Block  of  the  Test  Area.  Sub-blocks  #3  and  4. 
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Figure  48.  Comparison  by  Station  of  Control  and  Least  Squares 
Collocation  -  Predicted  Vertical  Deflections  at  the  South 
Block  of  the  Test  Area.  Sub-blocks  #5  and  6. 
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Figure  49.  Comparison  by  Station  of  Control  and  Least  Squares 
Collocation  -  Predicted  Vertical  Deflections  at  the  South 
Block  of  the  Test  Area.  Sub-blocks  47  and  8. 
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5.5.  Comparison  of  the  BJerhammar  and  Hardy  Predictora  to 
Least-Squares  Collocation 

In  the  sequel  a  comparison  between  the  Bjerhammar  Method  (BM), 
the  Hardy  Method  (HM)  and  Least  Squares  Collocation  (LSC)  will  be 
attempted.  The  CPU  time  comparison  is  based  on  times  required  on  the 
IBM  3081  computer.  In  the  event  that  only  Ag  observations  are 
available  the  results  from  the  White  Sands  tests  are  shown  in  Table  37. 


Table  37.  RMS  Differences  Between  Predicted  and  Control  from  the 
Three  Predictors.  Only  Ag  Observed. 


Block 

Sub-Block 

Aft 

[■gals) _ 

_ HV _ 

_ nil) _ 

CPU  Tine 

(sec) 

BM 

HM 

LSC 

BM 

HM 

LSC 

BM 

HM 

LSC 

NORTH 

3.08 

3.13 

2.84 

0.90 

4.93 

0.76 

1.00 

1L.36 

0.92 

BM 

SOOTH 

3.79 

3.88 

3.68 

0.90 

2.43 

0.80 

1.13 

4.77 

1.03 

NB:  512 

1 

4.27 

4.32 

4.23 

0.91 

2.86 

0.83 

0.81 

11.42 

0.79 

SB:  1440 

2 

2.54 

2.56 

2.59 

0.64 

6.17 

0.63 

0.99 

3.95 

0.91 

HM 

3 

2.84 

2.88 

2.28 

1.26 

2.32 

1.02 

0.70 

18.44 

0.75 

NB:  500 

4 

2.86 

2.98 

2.65 

0.81 

6.00 

0.61 

1.28 

9.52 

1.10 

SB:  1392 

5 

3.43 

3.73 

3.01 

1.27 

3.93 

1.00 

1.36 

7.24 

1.25 

LSC 

6 

6.22 

5.69 

6.25 

0.94 

1.49 

0.86 

1.28 

4.97 

1.16 

NB:  638 

7 

3.20 

3.57 

3.36 

0.58 

1.04 

0.57 

0.75 

4.18 

0.79 

SB:  1182 

8 

2.23 

2.63 

1.68 

0.74 

2.91 

0.67 

0.89 

3.12 

0.79 

Prom  Table  37  and  keeping  in  mind  that  both  the  observed  and  the 
control  Ag  have  a  standard  deviation  of  2  mgals  one  sees  that  all  three 
methods  can  predict  gravity  anomalies  within  3  to  4  mgals. 

Furthermore,  the  difference  in  the  quality  of  the  predictions 
introduced  by  each  method  never  exceeded  the  observation  error. 
Alao,  sub-blocks  that  performed  well  or  poorly  with  some  method 

behaved  similarly  with  all  methods.  For  instance,  sub-block  *6  yielded 
the  largest  RMS  difference  and  sub-block  #8  yielded  the  smallest  RMS 
difference  for  all  methods. 

The  picture  is  different  for  vertical  deflection  predictions.  From 
Table  37  one  can  see  immediately  that  the  HM  cannot  perform  to  a 
satisfactory  level.  On  the  other  hand  BM  and  LSC  performed  equally 
well  with  the  exception  of  sub-blocks  >3  and  5  at  which  LSC 
outperformed  BM  at  the  (  predictions  by  about  0125  which  is  not  very 

large  keeping  in  mind  that  the  standard  deviations  of  the  control 

deflections  are  013. 

In  case  that  both  Ag  and  ({,*))  are  utilized  as  observations  one 
gets  the  results  of  Table  38. 
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Table  38.  RMS  Differences  Between  Predicted  and  Control  from  the 
Three  Predictors.  Both  Ag  and  ({,»})  Observed. 
Downward  Continuation  on  a  (7'x7')  Grid  for  Both  BM 
and  HM. 


Block 

Sub-Block 

_ Ag(agals) _ 

_ tin _ 

_ tin _ 

mm\ 

BM 

LMJ 

LSC 

BM 

HM 

El 

mm 

i^«i 

NORTH 

3.56 

CE1 

2.93 

0.65 

0.81 

o 

CJl 

00 

0.65 

ESQ 

0.56 

BM 

SOUTH 

4.69 

5.34 

3.44 

0.62 

0.69 

0.64 

0.92 

1.08 

0.66 

NB:  892 

1 

4.07 

4.56 

4.08 

0.66 

0.69 

0.59 

0.54 

0.54 

0.48 

2 

3.16 

3.68 

2.40 

0.63 

0.62 

0.49 

0.75 

0.77 

0.60 

WMKMm 

3 

3.03 

3.85 

2.56 

0.83 

0.82 

0.44 

0,93 

0.31 

4 

4.34 

6.14 

EE3 

0.49 

0.44 

0.75 

0.79 

0.68 

5 

4.30 

4.77 

3.37 

0.57 

0.76 

0.89 

0.89 

1.18 

1.06 

LSC 

6 

4.06 

5.03 

4.61 

0.70 

0.75 

0.70 

0.80 

0.91 

0.59 

NB:  861 

7 

6.25 

BED 

3.70 

0.59 

0.69 

0.31 

1.03 

1.12 

0.37 

SB:  1469 

8 

2.90 

3.37 

1.82 

0.51 

0.55 

0.49 

1.04 

1.23 

0.59 

From  Table  38  one  sees  that  all  methods  can  predict  Ag  to  about  3 
to  5  mgals.  Also,  the  discrepancies  of  the  RMS  differences  among 
different  methods  are  always  smaller  than  the  standard  deviation  of  the 
Ag  values.  LSC  is  favored  over  both  BM  and  HM  in  terms  of  Ag 
predictions. 

As  far  as  vertical  deflection  predictions  are  concerned  all  methods 
can  predict  (  to  about  0!6  to  018  and  tj  to  about  016  to  1"0.  In  the 
majority  of  the  cases  all  methods  performed  equally  well  with  the 
exception  of  sub-blocks  #3  and  4  in  favor  of  BM  and  LSC,  #5  in  favor 
of  BM  and  #7  in  favor  of  LSC  for  |  and  sub-block  *3  in  favor  of  BM 
and  LSC,  *5  in  favor  of  BM  and  LSC  and  116,  7  and  8  in  favor  of  LSC 
for  i).  Overall,  similar  accuracy  was  obtained  by  all  three  methods. 

In  the  event  that  one  has  ((,q)  observations  only  one  gets  the 
results  of  Table  39. 

From  Table  39  one  sees  superiority  of  the  LSC  solution  in  the  Ag 
predictions,  which,  however  yields  rather  large  RMS  discrepancies.  In 
terms  of  vertical  deflection  predictions  one  can  observe  LSC  to  perform 
better  than  both  BM  and  HM  with  the  exception  of  sub-block  #5  for  tj. 
Therefore,  in  this  case  the  LSC  solution  iB  preferred. 
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Table  39.  RMS  Differences  Between  Predicted  and  Control  from  the 
Three  Predictors.  Only  Observed. 


Block 

Sub-Block 

Ag(mgals) 

_ iT) _ 

_ n<n _ 

■sgaa 

(gSI 

BM 

Kfl 

LSC 

BM 

MSI 

LSC 

LSC 

NORTH 

13.19 

16.48 

5.41 

1.25 

0.97 

(33 

1.06 

0.98 

0.86 

BM 

SOUTH 

9.19 

12.61 

7.36 

1.08 

0.91 

BBS 

0.99 

1.01 

0.86 

NB:  36 

1 

11.98 

11.68 

6.02 

1.06 

0.87 

0.65 

0.91 

1.16 

0.87 

SB:  36 

2 

8.15 

5.93 

3.75 

0.76 

0.74 

0.52 

0.74 

0.77 

0.66 

HM 

3 

16.00 

5.91 

ncin 

0.97 

0.73 

0.80 

0.88 

NB:  36 

4 

12.85 

5.46 

1.79 

1.20 

0.76 

1.50 

1.09 

1.02 

SB:  36 

5 

10.62 

6.49 

5.92 

1.58 

1.46 

1.31 

1.25 

1.23 

1.27 

LSC 

6 

7.34 

20.09 

3.49 

1.11 

0.58 

0.89 

0.93 

0.79 

NB: 130 

7 

10.20 

10.90 

10.70 

1.00 

0.97 

0.53 

1.35 

1.37 

0.64 

SB: 140 

8 

7.25 

12.46 

5.47 

0.75 

0.72 

0.60 

0.85 

0.86 

0.80 

A  comparison  of  the  three  methods  reveals  that  BM  and  LSC 
performed  equally  well  in  all  cases,  except  the  case  of  only  ((,*)) 
observations,  in  which  LSC  performed  better  than  BM.  The  HM, 
however,  yielded  peculiar  results.  For  example,  Ag  were  predicted  to 
about  3  to  4  mgals  from  gravity  data  only  (Table  37).  From  the  same 
solution,  (f,7j)  were  predicted  unacceptably.  From  only  (£,■>)) 
observations  (Table  39),  vertical  deflections  were  predicted  well, 
whereas  Ag  were  predicted  unacceptably.  These  types  of  results  from 
HM  come  as  no  surprise  in  light  of  the  comments  in  subsection  3.2.4. 

The  best  vertical  deflection  predictions  for  all  three  methods  were 
abtained  when  both  Ag  and  ((,»})  abservations  were  used  (Table  38). 
The  downward  continuation  for  both  the  BM  and  the  HM  was  performed 
on  to  a  (VxV)  on  the  geosphere.  In  Table  38,  only  the  RMS 
differences  of  control  minus  predicted  quantities  are  given.  The 
corresponding  average  differences  are  about  1  mgal  for  Ag  and  for 
(f,7j)  they  are  in  the  order  of  a  few  tenths  of  an  arcsecond. 
Furthermore,  the  predictions  obtained  from  ie  three  methods  agree 
very  well  as  seen  from  Figures  26  through  49.  The  RMS  prediction 
differences  between  any  two  methods  are  2  to  4  mgals  for  Ag  and  in 
the  sub-second  level  for  (f,7j).  Furthermore,  the  corresponding 
average  differences  are  less  than  0.7  mgals  for  Ag  and  for  ($,»j)  they 
are  less  than  two  tenths  of  an  arcsecond  in  absolute  value.  In  the 
North  Block,  the  best  agreement  is  observed  between  the  predicted 
both  Ag  and  (f,ij)  from  the  BM  and  LSC  and  the  worst  agreement  is 
observed  between  HM  and  LSC.  In  the  South  Block,  the  best 
agreement  is  observed  between  BM  and  HM  whereas  the  worst  one  is 
observed  between  the  HM  and  LSC.  Correlation  coefficients  between 
average  (f,»j)  differences  (control  minus  predicted)  among  methods  in 
each  sub-block  ranged  from  0.2  to  0.9.  The  average  correlation 
coefficient  was  0.7.  The  corresponding  correlation  coefficients  from  Ag 
predictions  ranged  from  0.82  to  0.99.  The  average  value  was  0.90. 
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Examination  of  the  differences  of  control  minus  predicted  vertical 
deflections  at  individual  stations  yields  interesting  results.  For 
example,  in  Kearsley  et  al.,  [1985,  p.  68]  the  f  component  of  the 
vertical  deflection  at  station  191  was  reported  as  a  suspected  error  in 
the  control  data.  This  appears  to  be  the  case  from  the  results  of  this 
investigation  also.  Furthermore,  stations  with  large  differences  from 
one  method,  yield  large  differences  with  all  three  methods  (e.g.  stations 
376,  383,  395,  404,  286  and  305  for  {  and  stations  199,  200,  202,  203,  128 
and  29  for  7)  to  name  only  a  few).  On  the  other  hand,  there  are  some 
stations  with  small  differences  from  one  method  and  large  differences 
from  another  (e.g.  stations  378,  320  and  172  for  (  and  349,  355  and  127 
for  ti). 

As  far  as  CPU  time  requirements  for  the  three  predictors,  Tables  37, 
38  and  39  indicate  that  there  is  no  method  that  consistently  required 
less  time  than  the  others.  It  is  worth  noting  that  75%  of  the  time 
estimates  for  BM  and  HM  is  needed  to  compute  an  optimal  value  for  the 
radius  of  the  geosphere. 

The  software  for  both  the  BM  and  the  HM  was  converted  to  work 
on  the  CRAY  X-MP/24  supercomputer.  This  conversion  was  almost 
effortless.  However,  it  will  take  a  moderate  effort  to  modify  the  LSC 
software  (GEOCOL)  to  work  on  the  supercomputer.  The  CPU  time 
requirements  for  both  the  BM  and  the  HM  on  the  CRAY  X-MP/24, 
including  an  optimal  r0  computation,  are  shown  in  Table  40. 


Table  40.  CPU  Time  Requirements  (in  seconds)  for  BM  and  HM  on 
the  CRAY  X-MP/24  Supercomputer. 


Type  of  Observable 

BM 

HM 

NB 

SB 

NB 

SB 

A* 

7.7 

22.6 

8.8 

23.7 

Ag  and  ($,ij) 

10.4 

24.8 

10.9 

28.8 

0.1 

0.1 

0.1 

0.1 

Comparison  of  Table  40  to  Tables  37,  38  and  39  reveals  improvement 
by  a  factor  of  at  least  10  and  as  much  as  90. 

A 

Lastly,  in  the  event  that  an  optimal  r0  is  required  for  the  BM  or 
the  HM  and  no  control  data  exist  in  an  area,  then  the  observations  can 
be  separated  in  two  groups.  The  first  group  should  be  regarded  as 
observations  and  the  second  one  should  be  regarded  as  control  data. 
These  two  groups  can  be  used  with  the  s3-method  to  compute  an 
optimal  geosphere  radius.  Finally,  the  entire  data  set  should  bo  used 
as  observations  to  perform  the  solution. 
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5.6  Comparison  With  the  Pour  Methods  Tested  with  the  New  Mexico  Test 
Data 

In  Kearaley  at  al.  [1985]  four  methods  to  predict  deflections  of  the 
vertical  from  gravity  anomaly  data  were  tested  and  intercompared. 
These  methods  were  the  Fast  Fourier  Transform  (FFT),  the  Combined 
Collocation-Integration  (CINT),  the  Numerical  Integration  (RINT)  and  the 
Terrain  Effect  Integration  and  Collocation  (TEIC)  method.  In  Tables  4.5 
and  4.6  of  [ibidi  pp.  93-94]  they  reported  RMS  discrepancies,  control 
minus  predicted  vertical  deflections,  predicted  from  gravity  data  in  the 
order  of  1"  when  height  data  are  used. 

In  the  sequel  we  will  compare  the  results  of  the  four  methods  in 
[Kearaley  et  al.,  1985]  with  the  results  of  this  investigation.  The 
comparison  will  be  based  on  Tables  4.5  and  4.6  of  [ibid,  pp.  93-94]  and 
on  Table  37  of  this  work.  However,  the  results  of  our  LSC  will  be 
used  rather  than  the  ones  of  TEIC  because  the  results  of  our  LSC 
solution  are  slightly  better.  A  summary  of  these  results  appear  in 
Table  41.  In  Table  41,  the  column  designated  AG  shows  the  RMS  values 
of  the  vertical  deflections  by  sub-block. 


Table  41.  Comparison  of  the  Bjerhammar  and  Hardy  Predictors 
with  the  Four  Methods  Tested  at  New  Mexico.  Only  Ag 
Observed. 


AG 

BM 

m 

LSC 

FFT 

CINT 

RINT  1 

BLOCK 

n 

n 

09 

ism 

09 

in 

09 

09 

09 

09 

09 

09 

09 

09 

4.3 

7.9 

w 

IK] 

cun 

11.4 

0.8 

HR 

i.i 

IKI 

IKI 

Bid 

l.l 

2.0 

2.8 

5.6 

EE 

0.8 

2.9 

11.4 

0.8 

0.8 

1.2 

1.4 

0.8 

1.4 

0.7 

1.4 

2 

4.1 

9.1 

BE 

•s 

6.2 

Kj 

EE 

1.0 

1.4 

0.9 

1.6 

0.7 

2.2 

3 

5.1 

6.6 

1.3 

8 

2.3 

18.4 

m 

1.2 

1.3 

1.4 

1.5 

1.8 

1.6 

4 

4.7 

8.7 

0.8 

m 

I3M 

•M; 

i.i 

0.8 

1.3 

ran 

1.8 

0.9 

2.5 

South 

m 

7.1 

□E 

1.1 

2.4 

4.8 

& 

jrn 

EE 

1.2 

0.8 

1.2 

0.7 

1.6 

5 

7.9 

1.3 

1.4 

3.9 

7.2 

£ 

m 

1.2 

0.7 

1.3 

0.9 

1.6 

6 

2.8 

7.7 

1.3 

1.5 

E2y 

0.8 

l.l 

0.9 

1.1 

0.6 

1.7 

7 

6.3 

0.6 

PJd 

be 

4.2 

0.6 

0.8 

m«i 

1.4 

0.6 

1.2 

8 

Rfl 

MM 

lUJ 

YW 

3.1 

F9i 

m 

0.8 

1.1 

0.8 

1.2 

0.7 

W3 

From  Table  41  one  can  see  that  the  methods  that  performed  best 
were  LSC  and  BM.  Furthermore,  LSC  performed  slightly  better  than 
BM  (by  about  (711  on  the  average,  which  is  below  the  ((,y)  noise 
level). 

However,  the  most  important  conclusion  drawn  from  Table  41  is 
that,  when  reference  field  and  RTM  effects  are  removed  from  the  data 
and  restored  at  the  predictions,  there  are  at  least  five  methods  that 
can  predict  vertical  deflections  to  the  sub-second  level  from  gravity 
data. 
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CHAPTER  VI 


SUMMARY,  CONCLUSIONS,  RECOMMENDATIONS 


Two  deterministic  methods  for  gravity  field  approximation  have 
been  investigated.  The  first  one  was  the  Bjerhammar  Dirac  Impulse 
method  and  the  second  one  was  the  Hardy's  biharmonic  potential 
method. 

Bjerhammar  defined  the  discrete  geodetic  boundary  value  problem 
as  the  one  at  which  observations  are  given  at  discrete  points  and  it  is 
required  to  find  a  gravity  field  Buch  that  all  observations  are 
satisfied.  The  solution  is  constructed  such  that  the  disturbing 
potential  is  harmonic  outside  a  sphere  fully  internal  to  the  Earth  and 
regular  at  infinity.  The  Dirac  Impulses  that  generate  the  disturbing 
potential  are  computed  by  a  downward  continuation  process  and  they 
are  used  to  perform  predictions  in  an  upward  continuation  scheme. 

Hardy's  work  was  initiated  by  the  fact  that  the  integral 
representation  of  the  disturbing  potential  is  singular  at  points  that 
induce  potentiaL  Based  on  the  non-uniqueness  of  the  solution  of  the 
inverse  problem  of  potential  theory  one  can  assume  that  the  density 
anomaly  function  of  the  Earth  together  with  its  normal  derivative  is 
zero  at  the  boundary.  An  integral  representation  of  the  disturbing 
potential  can  be  derived  which  is  non-singular  at  points  that  induce 
potential  and  which  satisfies  the  biharmonic  equation.  An 
approximation  to  the  fundamental  integral  is  also  derived. 
Operationally,  the  biharmonic  sources  are  computed  as  the  solution  of  a 
linear  system  and  then  they  are  used  to  perform  predictions. 

Both  of  the  aforementioned  methods  can  use  any  linear  functionals 
of  the  disturbing  potential  as  observations  and/or  quantities  to  be 
predicted. 

Tests  were  performed  for  both  methods  with  the  White  Sands  Test 
Data.  The  predictions  were  compared  to  independently  observed 
gravity  anomalies  and  vertical  deflections  that  served  as  control  data. 
Reference  field  and  residual  terrain  model  effects  were  removed  from 
the  observations  and  were  restored  at  the  predicted  values  before  any 
comparison  to  the  control  data  was  done. 

A  factor  that  influences  the  quality  of  the  results  with  both 
predictors  is  the  radius  of  the  internal  sphere.  Two  approaches  to 
compute  it  from  gravity  data  failed  for  both  methods.  However,  a 
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technique  for  optimal  radius  computation  that  yielded  satisfactory 
results  is  to  consider  some  measure  of  the  quality  of  the  predictions  a 
second  order  polynomial  in  the  radius. 

The  Bjerhammar  method  (  performing  the  downward  continuation  on 
the  nadir  points  of  the  observations  and  with  only  gravity  data 
resulted  in  RMS  differences  of  control  minus  predicted  values  in  the 
order  of  3  to  4  mgals  Ag,  (719  for  (  and  110  tor  rj.  When  vertical 
deflection  observations  were  introduced  the  RMS  discrepancies  became 

3.5  to  4.5  mgals  for  Ag  and  018  for  (  and  q.  When  the  gravity 
observations  were  completely  removed,  the  RMS  differences  became 
larger  than  10  mgals  for  Ag,  and  about  1"  for  both  (  and  ij.  The 
aforementioned  results  pertain  to  both  the  Asymmetric  and  the 
Symmetric  Kernel  approach.  In  fact,  the  only  difference  between  the 
AK  and  the  SK  is  that  the  optimal  radii  associated  with  the  AK  are 
usually  smaller  than  the  ones  with  the  SK  (see  Table  17). 

When  the  downward  continuation  is  performed  onto  a  grid  on  the 
geosphere,  the  Bjerhammar  method  predicted  Ag  to  3  to  4  mgals,  {  to 
about  019  and  ij  to  about  I'll  from  gravity  observations  alone.  When 
vertical  deflection  observations  were  introduced,  the  RMS  differences  of 
control  minus  predicted  values  was  the  same  (about  3  to  4  mgals)  for 
Ag,  whereas  it  became  about  017  for  (  and  018  for  ij.  On  the  other 
hand,  the  predictions  from  ((,»})  data  alone  were  unacceptable  both  for 
Ag  and  (£,»,). 

For  every  variation  of  the  Bjerhammar  predictor,  standard 
deviations  of  the  predictions  were  computed  according  to  (2-57).  These 
standard  deviations  cannot  be  considered  a  safe  indicator  of  the 
quality  of  the  predictions.  The  reason  for  this  is  that  there  were 
many  well  predicted  quantities  with  large  standard  deviations  and 
there  were  many  poorly  predicted  quantities  with  small  standard 
deviations. 

The  Hardy  method,  when  the  biharmonic  sources  were  located  at 
the  nadir  points  of  the  observations,  gave  RMS  discrepancies  of  control 
minus  predicted  values  on  the  order  of  3  to  4  mgals  for  Ag  whereas 
the  vertical  deflection  predictions  were  very  poor.  When  vertical 
deflection  observations  were  introduced,  the  RMS  differences  were 
larger  than  10  mgals  for  Ag  and  larger  than  215  for  (  and  rj.  When  the 
gravity  data  were  completely  removed,  the  RMS  differences  were 
degraded  further  tor  Ag  whereas  they  became  smaller  than  1"  for  both 
{  and  ij. 

When  the  downward  continuation  is  performed  onto  a  grid  on  the 
geosphere,  the  Hardy  method  with  only  gravity  data,  yielded  RMS 
discrepancies  in  the  order  of  3  to  4  mgals  for  Ag  whereas  the  vertical 
deflections  were  worse  than  217  for  all  grid  sizes.  The  introduction  of 
vertical  deflection  observations  degraded  the  Ag  predictions  to  4.5  to 

5.5  mgals  whereas  it  upgraded  the  ((,7j)  predictions  to  the  1"  level. 
The  complete  removal  of  Ag  data  rendered  both  the  Ag  and  the  (f,n) 
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predictions  unacceptable. 

For  every  variation  of  the  Hardy  predictor,  the  standard  deviations 
a  of  the  predictions  did  not  prove  to  be  indicative  of  the  quality  of 
the  resultss,  since  in  many  cases  large  a' a  were  associated  with  well 
predicted  quantities  and  vice  versa. 

Comparison  of  the  two  predictors  (BM  and  HM)  with  Least  Squares 
Collocation  (LSC)  indicates  that  BM  and  LSC  yield  comparable  results  in 
all  cases  with  the  exception  of  the  case  where  only  (f  ,»j)  observations 
are  utilized  in  which  case  LSC  performed  better  than  BM.  On  the 
other  hand,  HM  performed  well  when  it  predicted  Ag  from  Ag  data  or 
from  {{,?))  data  and  the  downward  continuation  was  performed  at 
the  nadir  points  of  the  observations. 

The  most  important  overall  result  of  this  work  is  that  when 
reference  field  and  RTM  effects  are  taken  into  account,  there  are  at 
least  five  methods  that  can  predict  ((,ij)  from  Ag  to  the  sub-second 
level,  even  in  mountainous  areas.  Furthermore,  the  improvement  of  the 
predictions  should  not  be  anticipated  from  a  theoretical  breakthrough 
but  from  data- type  and  coverage  improvement. 

As  far  as  future  investigations  are  concerned  it  is  recommended 
that  undulations  and/or  gravity  gradients  be  predicted  from  various 
data  types  with  both  predictors.  Particularly  for  the  Hardy  method  it 
is  suggested  that  a  low  degree  and  order  spherical  biharmonic 
expansion  (e.g.  6  to  10)  be  computed  from  the  formulae  given  in 
Appendix  A.5  using  10'xlO*  or  5*x5*  global  data  and  then  be  tested  as 
to  its  reliability. 


APPENDIX  A 


DERIVATIONS 


1.  Show  that  if  M,  = 

with  bf  =  t  -  3dt  +  -  5tacosu  -  3tacosuinu, 

t*-l 

then  M,  =  1  +  +  3tcosu 

Proof: 

Recalling  equations  (2-21)  and  (2-22)  one  has 

3d  t-coscj 
3t  ~  d 

and 

3u  1  f  t-CQScj 

at  '  2  l  d 

Therefore 

Mj  =  ~ ^  =  2  -  6d  +  j  -  lOtcosu  -  6tcosuinu  - 

~(l  _  (td)  +  -  lOtcosu  -  6tcosu<nu  -  3t2cosu  |^(<nu) 


-  COSCJ 


). 


=  2  -  6d  +  j  -  lOtcosu  -  6tcosuinu  -  |l  -  3t  -  -C°S{J  -  3d  + 

•  .  t-cosu 

^  *  d 

+  2  -  -  lOtcosu  -  6tcosuinu  - 

-  3t2cosu  i  I(1Z£22"  _  COSMjJ  = 

=  2  -  6d  +  4  -  lOtcosu  -  6tcoswinu  -  1  +  “  — — PSu  +  3d  - 

a  ad 
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+  _  2tcosu  +  iotcos^Btco^inu  +  ^.(y^t-coso)  _  3t?|osau  = 

d9  d9  2ud  2u 

_  1  -  3d  +  —  +  _  3tcosu  +  2ta  _  2tcosu  +  3taco8u(t-co8u)  _ 

d  d  d  ds  d9  2ud 

3tacosau 

2u 

,  «,  .  2  .  3t*  6tcosu  .  2ta  2tcosu  .  3tcos<u  . 

-  1  "  3d  +  d  +  “d - d -  +  "d*  ‘  +  - d~  + 

3tacosu(t-cosu)  _  3tacosau 
2ud  2u 

=  jy  [d9  -  3d*  +  2d*  +  3t*da  “  6tdacosu  +  2ta  -  210086,]  + 

3tCOS6>  r.  t a- tCQ86r-tdCQ86>l 

d  l1  2u  J  = 

=  -ji  [ds  -  3d4  -  da  +  3da(l  +  ta  -  2tcos6>)  +  2ta  -  2tcosu]  + 

+  [2u  +  ta  -  tcosu  -  tdcosu]  = 

=  jy  (d9  -  3d4  -  da  +  3d4  +  2ta  -  2tco86»)  + 

+  (1  -  tcosu  +  d  +  ta  -  tcosu  -  tdcosu)  = 

=  jy  (d9  -  da  +  2ta  -  2tcosu)  +  3t^“  (da  +  d  -  tdcosu)  = 

=  jy  (d*  -  1  -  t*  +  2tcosu  +  2ta  -  2 tcosu)  + 

+  3t^"  d  (d  +  1  -  tcosu)  S 

=  jjy  (d*  +  ta  -  1)  +  2u  =  1  +  “pp  +  3tcosu  q.e.d. 

2.  Show  that  A<  =  [l  “  3d  +  4  -  5tcosu  -  3tcosu4nu]  = 
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■  ‘"‘""l8  -  b  -  * 3,nu) 


Proof: 

Recalling  equations  (2-16)  and  (2-18)  one  can  show  that 

Therefore: 

A|  =  [l  -  3d  +  g  ~  5tcosw  -  3tcoacjinuj  = 

=  -3  +2(-l)d-a  -5t(-sinu)  -3t(-sin«)4nu  -3tcosu  ~  = 


=  -3 


tsinu  2  tsinu 
d  ~  da  d 


+  5tsinu+3tsinuinu  -  3tc°s-  +  j)  = 


*  2  -  fj  ♦  5  ♦  3100  -  -  2i^i]  = 

=  tsinu ^8  ~  js  +  3Anu  -  (2ud  +  2u  +  tdcosu  +  tcosu)]  = 


=  tsinu 


in«[8  - 


+  34nu  - 


2^|  ((1  -  tcosu  +  d)d+(l  -  tcosu  +  d)  + 


+  tdcosu  +  tcosu) ]  = 

12  3 

8  -  +  3<nu  -  gj—  (d  -  dtcoscj  +  d*  +  1  -  tcoscj  + 

+  d  +  tdcosu  +  tcosu)]  a 

=  tsinu [8  -  —■  +  31nu  ~  ^  (da  +  2d  +  1)]  = 

=  tsinu ^8  -  +  3inu)  q.e.d. 

3.  Show  that  if  gff  =  -  3t3cosu  -  ta, 

with  da  =  1  +  ta  -  2tcosu  and  t  = 

r  ’ 

then  fjU  =  “  [3?  (1  -  2ta)  -  -  9tcosu  -  2] 
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Proof: 

Firstly:  t  =  =>  (t")  =  nt""‘  =  nt«-»  i  =  l£l 

r  *ro  *r0  r  r0 

Secondly:  [<l  +  t*  -  2tcos«)*]  = 

=  <2t-2°»">  £  ■ 

Hence: 


h>.  Igggg:  *  *’^1-  ~  *•»■*• 

ar«  d* 


d* 

-  3co~  . 

*r„  3r0 


a _ 


*  a*  (a-**)  “  *  **  (-  *£)]  -  ^  - 


■  3coso 


3ts  2ta 


a  (1  ~  2ta)  -  ^>-^1~-ta)(t~coa<J)  _  cose  _  2ta 


r0d* 


r0 


,  or 


Mi ,  n  ra_  „  3tu-ta)(t-cos«)  _  _i 

»r0  r0  Id*  (1  Zt  '  d*-'  1  ~  9tcos«  -  2j  q.e.d. 

4.  Show  that  if  f  flj  ]  =  SiSSl  fcos*|  , 

'  ’  7  l»inccl  * 

i .  st*  -  %  -  Kgu:  *  3f«M, 

u  =  5  (1  ~  tcosw  +  d)  and  I  and  t  as  in  A. 3,  then, 


ai* 


£1  “S*  ♦ 


+  ltlTC°»W  (2  jd+1)*  _J_)  _  tcose  f(d+l)a  .  1 1  f cose) 

d  U*  +  4ua  +  2udJ  l  2Sd  +  MJ  (sinaJ 


1 18 


Proof : 


Firstly:  £  ■  |  [-  cm.  *  £]  = 


1] 

_  t 

t-coso  1 

2r0 

JJ -  COScjj 

Then: 


d»  3(da) 

(a)  (£1)  =  - 2Eh _ ! _ ilfl.  _  3t®  t*-3d*  t(t-cos«) 

,r*  ld’j  d*  -^sr  “7^ — ’ 

_i_  fill  _  3ti_  _  3t4 (t-coaa) 

3r0  Id* J  r„ds  r0ds 

(b)  fc  fts(d+l)*]  =  (d+1)*  iiiii  +  t *  [(d+1)*]  = 


or 


=  (d+1) *  311  +  ts -2(d+l) 


ad 


=  [3(d+D  +  S&rP—”)] 

4 8t ♦  •  ft  ■ ♦ . itopd. 


thus 


a _  rt*(d+p*i  .  ud  T77  (t»(d+l)*]  -  t*(d+l)* 

ar0  l  ud  J  Jaja  1 

1  i^51  13<W)  +  - 

- [&  P3**  -  —)  *  *“&■*] 

(')  £  [tMm.]  =  jnu  iUil  *  t>  (l„„)  = 

=  211  i~.,  ^  li  t  ft-coau  1 

r.  ,nu  +  “  2^  l~ - <*>•«] . 


(ud) 


therefore: 
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2t4(d»lHt-cos«)  t4(d+l)  * 
r©uda  2r0uada 


(t  -  cosu  -  dcoscj)  - 


_  t4(d+l)a(t-coep)l  +  9t9  ,  +  3t4(t-cosw)  _  3t4coa<a  _ 

r*ud*  J  r0  2r0ud  2r0u 

=  —  (at9  -  +  3t*<nu)  +  ^^5§Sggl  - 

r©  l  d9  2ud  /  r0d* 

3t4(d»lHt-coat»)  .  3t4(d-*-l)a(’ t-coatr-dcoao) 
r«uda  4r0u*d* 

+  3t4(d+l)a(t-coa«>)  +  3t4(t-co»6>)  _  3t*cos<j 
2r«ud9  2r0ud  ~  2r0u 


3t4cos« 

r(d+na  .i 

4ua  J 

2ur0 

L  2Sd  +  1J 

(^r) +  l\'  thus: 


+ 


\tL 

»r« 

i£h 

*r© 


3t9»inw 

Tor 


( - 1.  -  ^  ♦ 


+  ^(Jtssssl  f2_  +  (dfl)a  JJ  _  tcoaw  f(ch-l)a  .llfccs*) 

da  Id9  +  4ua  +25di  l  2ud  NilainJ 


■ino(J 

q.e.d. 


5.  The  biharaonic  equation  is  AaV  =  0.  Find  ita  solutions. 

Solution: 

let  the  Cartesian  rectangular  coordinates  x,  y,  z  be  expressed  as: 

x  =  x(q»,  qa,  q,),  y  =  y<q»,  q*,  q,)t  z  =  z(q,,  qa,  qs) 

such  that  x,  y,  z  are  continuously  differentiable  functions  and  also 
solvable  for  q,,  q2,  q,,  i.e.,  the  Jacobian  of  the  transformation 
does  not  vanish. 


For  orthogonal  coordinate  systeas  [4^-  ill-  +  - 

*  7  Uqi  *qj  *q<  *qj  *q<  *qj 

0,  i  *  j] ,  the  Laplacian  ia  [Kellogg,  1929;  p.  183] ,  [Heiskanen 


where 


:,  1967;  p. 

19]: 

1  fJ 

fh.h. 

JV  1  . 

J 

hihfhs l Jq, 

l"b^ 

Jqj  + 

»q« 

/(*-]*  + 

|ii_r 

+  fi£-] 

n* 

UjqJ 

IjqJ 

ljq,J 

i  j 

and  hthshs  =  I f — rl  =  IJJ  and  J  ia  the  Jacobian  of  the 

|a(q».qi,qs)l 

transformation. 

In  the  usual  spherical  coordinates  which  satisfy: 


x  =  rainflcoaX 
y  =  rsinOsinA 
z  =  rcosfl 


one  has 

Jx  _ 
Jr 

sinOcoaX; 

Jx 

99  ~ 

rcosflcosX; 

Jx 

JX 

II 

sindsinX; 

II 

rcosOsinX; 

II 

Jz 

Jr  " 

cos0; 

Jz 

99  ' 

-  rsin0; 

Jz 

JX  " 

h,  =  1;  h,  =  r;  hs  =  rsin0. 

With  these  values  for  q<,  h«,  i  =  1,  2,  3,  AV  becoaes  after 
differentiation  [Heiskanen  and  Moritz,  1967;  p.  19]: 

lv  Jav  ,  2  jv  ,  l  j»v  cotg  jv  l  j»v 

Jr*  r  Jr  r*  J0*  r*  99  rJsin *9  JAa 


**'  -  “**>  *  pfej  [fc  [H—  *  IS  (— 

*  i*  fcA  « ■  pin  Gi  ♦  3  ♦  ffl 
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Computing  A,  B,  f  yields: 


a  isv  iu  jin 

A  =  rasin©  —  (AV)  =  rasin®  jp?  -  2a in©  —  +  2rsin®  ypj 


2sin 0  3aV  .  .  a*V 

T”  Jo*  +  8infl  a^e* 


geos©  JV±  .  a2V 
T"  a®  +  cos*  arae 


a*v  +  l 


a*v 


rsin®  a\2  sin®  3r3Xa 


B  =  sin®  jg  (AV) 


sin® 


a2v 

aaar* 


+  2sin®  a2V  +  sin®  3*V 


aear 


r2sin2®  a© 


av  +  cos®  a2v 


2cos®  aaV 


a®* 


a*v 


r2  a®2  r2sinae  ax2  rasin®  aaax2 


r5  J.L  f  AV)  =  -i-  ili-  +  -?-■  i!3L  +  _J _ 2i!L 

sin®  ax  '  '  sin®  axar2  rsin®  axar  r2sin®  3X3©2 

+  -^£2 a«  +  -j-L—  1IV  therefore: 

rasina®  3X38  rasin*a  ax*  * 


AaV 


=  fiiv  +  l  a^  + 


1 


a4v  +  2 


a4v 


a4v 


lar4  r4  aa4  r4sin4a  ax4  ra  ar2aea  r2sin2a  araax2 
2 _  a4v  4  ayy  +  2cota  a2v  2cota  aav 


r4sin2a  aeaax2  r  ar*  ra  araae  r4sin2a  aaax2 

2cote  a*v  _  cot 2 a  aav  4  a»v  coteq+2sin2®)  av] 

r4  aa*  r4  aa*  r4sin4a  ax2  r4sinaa  aai 

Assusing  V(r,  0,  X)  =  f(r)'Y(8,  X)  with  Y(a,  X)  =  g(©)-h(X)  we  get 

A.y  =  f(«)Y  +  L.  ill  t  +  1  f  W  +  2_  f. .  i!Y  + 

r  Y  r4  aa4  f  r4sin4a  f  9\*  r2  r  aa2 


2  a*Y  2 

r2sinaa  f  ax2  +  r4sina9  *  aa^ax 


-  4.  4  ,(,)  v  . 

1  aa2ax2  r  1  1 


2cote  ...  ay 
ra  1  aa  “  r4sina® 


2cota  .  a*Y 
1  aaax2 


2cota  a*y  .  cot 2 a  .  a^Y  4  r  aaY 
r4  aa*  T  "  r4  1  aa2  r4sin4a  ax2 


cot8(l+2sina8)  ,  3Y 
r4sin2a  1  aa 
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where  the  derivatives  of  f  are  with  respect  to  r. 
Now  A*V  =  0  <=>  ~  A*V  =  0  <=> 


r4f(4)  .  1  J4Y  ,  1  a4Y  .  2raf 1  ,  2f"r*  1  a*Y 

f  +  y  ae4  +  Ysin4e  ax4  +  f  y  ae*  f  Ysin*a  axa 

2  a4Y  4r*f<s)  .  2r*f"  cot®  9Y  2cot8  aaY  , 

Ysinaa  ae*ax*  +  f  +  f  y  ae  "  Ysinaa  aeaxa 

2cota  a*Y  cotaa  aaY  ,  4  aa y  cote(i+2»in*e)  ay  _  0  m 

y  ae3  y  aaa  Ysin4e  axa  Ysinaa  ae 


...  aaY  ^  4  aaY  .  cote(H2sinaa)  9Y 

001  9  aaa  sin4a  axa  sin*e  aa 

r4f (*)  4r3f<3> 

f  f 


On  the  other  hand: 


az  _  a_ 
aa  "  aa 


aaY 

aea 


+  cotS 


£Y 

aa 


+ 


l  aaYi 

sinaa  ax*J 


asY  ,  .  _  aaY  l  ay  1  a3Y 

“  aa3  aa*  “  sinaa  aa  sinaa  aaax* 


2cosa  aaY 
sin3®  ax*’ 


and: 


aaz  a  ra3Y  .  aaY  l  ay  l  a3Y  2cos8  a*Yl  _ 

aaa  =  ae  lae3  cot®  aaa  *  sinae  ae  sinaa  aaax*  sin3e  ax*J 


a4y  _  l  aay  . ,  a3Y  .  2cosS  ay  _  l  aaY 
=  ae4  sinae  aea  cotw  ae3  sin3a  ae  sinae  ae* 


2  cos  a  a*Y  1  a4Y  2n+2cos*e)  a*Y  2  cose  a3Y 
sin3a  aaax*  sinae  aeaaxa  sin4a  ax*  “  sin3e  aeaxa 
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az  i»y 

IX  "  1X10* 


+  cotfi 


a»Y 

1X10 


l  i*y 

sin*0  IXs 


and 


l*Z 

IX* 


i4y  a*Y  l  i4y 

I0*IX*  +  cot®  I0IX*  +  ain*0  IX4 


hence: 


A,Z  +  2Z  =  0  + 


«  l  afz 

0010  10  ain*0  IX* 


l4Y  1  l*Y  -  afY  2coa0  IY  1  l»Y 
I04  "  ain*0  10*  10s  ain*0  10  "  ain*0  10* 

2co#0  l*Y  1  l4Y  2(l+2coa*0)  l»Y 
"  ain*0  I0IX*  ain*0  I0*IX*  ain40  IX*  ' 


2coa0  l»Y 
•in*0  I0IX* 


+  cot0 


a*Y 

10* 


+ 


cot0 


l*Y 

10* 


__L_  il_  + 

sin*0  10 


.  1  l»Y  2coe0  l*Y)  1  f  l4Y  l»Y 

sin*0  I0IX*  “  ain30  IX*i  sin*0  ll0*lX*  cot®  I0IX* 


+  _i—  i!2) 

sin*0  IX4J 


l*Y  IY 

2  a®?  +  2cot«  + 


2  l»Y 
sin*0  IX* 


34Y  ,  1  l4Y  2  l4Y  „  . _  IfY  2cot0  l*Y 

~  I04  ain40  IX4  ain*0  I0*IX*  °  10*  “  ain*0  I0IX* 


+ 


ain*0 


aliil  ■ f2+4coa*0 
*J|0*  l  ain40 


2cos»0  2  ]  1^ 
ain40  sin*0J  IX 


fcoaO  2coa01  IY  _ 

lain*0  ain0  i  10 


_  I^Y  1  l4Y  2  l4Y  ,  to  i!l  2cot0  l*Y 
'  I04  ain40  IX4  ain*0  I0»IX*  ™  10*  "  ain*0  I0IX* 


-  cot*0 


l*Y 

10* 


4  a*Y 
ain40  IX* 


cot0(H-2ain*0)  IY 
ain*0  10 


=  L 


therefore  AtZ  +  2Z  =  L 
Now  (1)  becoaes: 


2r*f '  * 
f 


0  <=> 


1  2r*f ' *  7 

<=>  M  +  ±  (4lZ  +  2Z)  +  £  =  0  <=> 

<=>  M  +  i  (AfY  +  2A| Y)  +  2r*y  ^  =  0 


(2) 
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If  one  assures  A>Y  -  cY  (3) 

with  c  a  constant  one  gets 

A?Y  =  A,(AtY)  =  At(cY)  =  cAtY  (4) 


Substituting  (4)  in  (2)  one  obtains 

1  O-lf' '  *  V 

M  +  ^  (cAiY  +  2AjY)  +  =y - =  0  which  using  (3)  becomes : 


1  2r*f ' 

N  +  ±  (c*Y  +  2cY)  + 


c  =  0  <=> 


<=>  ^  +  4r^<0  +  2rX_  c  t  +  2c  .  „ 


(5) 


Selecting  the  constant  c  =  -  n(n+l),  equation  (5)  yields 

HI  +  4  iff  _  2n(nfl)  a^f  +  n(n-l)(n+l)(«n-2)  f=Q  (6) 

ar4  r  ar*  r*  »ra  r* 

The  solutions  of  (6)  are  r",  r”(B+1),  r"+*  and  r~(n~l)  as  can  be 
verified  by  substitution.  On  the  other  hand  (3)  becomes: 

jiv  jy  1  ]>v 

W  *  cote  Jg  +  +  n(n+l)Y  =  0 

the  solutions  of  which  are  [Heiskanen  and  Moritz,  1967;  p.  21]  the 
surface  spherical  harmonics 

Pna(co*6)cosaX  and  Pna(cos0)sinaX 

and  Pna(cos0)  are  the  Legendre’s  functions  [ibid,  p.  21].  Therefore, 
the  biharmonic  functions  (solutions  of  A*V  =  0)  can  be  represented  as: 

V,(r,0,X)  =  Z^r"  So[(a„w  +  r*cna)cosaX  + 

+  (Km  +  r#dna)sinnX]Pna(cos0) 

•  In 

V.(r,0,X)  =  z  -JTT  I  [(a«»  +  r*cna)cos«X  + 

n=o  r  «=o 

+  (b»w  +  r*dna)siis«X]Poa(cos0) , 
where  a„a,  bna,  cna,  and  d^a  are  arbitrary  constants. 
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